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Fast, reliable orbital evolutions of compact objects around massive black holes will be needed as 
input for gravitational wave search algorithms in the data stream generated by the planned Laser 
Interferometer Space Antenna (LISA). Currently, the state of the art is a time-domain code by 
[Phys. Rev. D81, 084021, (2010)] that computes the gravitational self-force on a point-particle in 
an eccentric orbit around a Schwarzschild black hole. Existing time-domain codes take up to a few 
days to compute just one point in parameter space. In a series of articles, we advocate the use 
of a frequency-domain approach to the problem of gravitational self-force (GSF) with the ultimate 
goal of orbital evolution in mind. Here, we compute the GSF for a particle in a circular orbit in 
■ Schwarzschild spacetime. We solve the linearized Einstein equations for the metric perturbation in 

Lorenz gauge. Our frequency-domain code reproduces the time-domain results for the GSF up to 
, ~ 1000 times faster for small orbital radii. In forthcoming companion papers, we will generalize 

(N- our frequency-domain computations of the GSF to include bound (eccentric) orbits in Schwarzschild 

' spacetimes, where we will employ the method of extended homogeneous solutions [Phys. Rev. D 78, 

, i 084021 (2008)]. We wiU eventually extend our methods to attempt a frequency-domain computation 

of the GSF in Kerr spacetime. 



I. INTRODUCTION 



j^jj With the start of the upgrades to second generation ground based gravitational wave detectors and the approval 
of the LISA Pathfinder mission the age of gravitational wave (GW) astronomy has begun. One promising source 
of gravitational radiation is the so-called extreme mass ratio inspirals (EMRIs) where a compact object (a black hole 
^ . or a neutron star) of a few solar masses slowly spirals in toward a massive black hole (MBH). The compact object 

■ (CO) interacts with its own gravitational field, which causes it to move on a path perturbed from the geodesic of the 
\^ , background spacetime. Along this 'forced' trajectory, the object radiates gravitationally losing energy and angular 
00 ■ momentum. For CO to MBH mass ratios of ~ 10~^ — 10~^, the frequency of the gravitational waves emitted during 
I/"} ] the last few years of inspiral (up to the final plunge) will be a few mHz, which will fall right in the middle of LISA's 

■ frequency band 0. Analysis of the waveforms emanatingfrom these inspirals will provide us with an unprecedented 
way of mapping spacetime around the central objects [7|, which are presumed to be Kerr black holes. A typical 
LISA bandwidth EMRI will be a ~ 1.5Mq neutron star/black hole inspiraling onto a ^ lO^M© MBH. In its last year 
before the plunge, the compact object will spiral in from a distance of ^ lOGM/c^ to the innermost stable circular 
orbit (6GM/c^ for Schwarzschild black hole) executing ^ 5 x 10^ orbits and sweeping the GW frequency band from 
~ 2 mHz to ~ 5 mHz Q. Such sources will be detectable by LISA for years, but the amplitude of the resulting 
gravitational wave strain will be smaller than the noise in the instrument Q. However, matched- filtering the signal 
over an extended period of time few years) will bump the signal-to-noise ratio as high as 100 for the nearest sources 
[loj . To be able to use matched- filtering, very accurate gravitational wave templates will be required as input for the 
cross-correlation. This will call for very accurate simulations of these inspirals over their LISA bandwidth lifetimes. 
The most challenging part in obtaining reliable simulations will be keeping track of the orbital phase as over the 
course of the inspiral the accumulated phase error should not exceed a few radians out of a total of 0(10^) — 0(10^) 
radians. This will put quite a stringent limit on the error tolerance of orbital evolution models. 

This is where the gravitational self-force comes in. In the test mass — 0) case, the CO follows a geodesic of 
the background spacetime. However, for a small, but finite mass the CO (modelled as a point particle sourced by 
a Dirac delta function) interacts with its own gravitational field, which scatters off the curvature of the background 
spacetime. This interaction can be interpreted as perturbing the particle's path off the background geodesic. In other 
words, the particle now accelerates, thus feels a net force due to this back-reaction. This is what has become known 
as the gravitational self- force (GSF). 

The study of radiation reaction began not with the GSF but with electromagnetic self- force (SF). This problem 



X: 



was first successfully worked out by DeWitt & Brehme [11|. Later, the solution to the gravitational problem was 
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formulated by Mino, Sasaki & Tanaka [12| and independently by Quinn & Wald [13| in terms of "forced geodesies" 
where the compact object feels a net force and is pushed off the geodesic of the unperturbed background spacetime. 
This approach is generally known as the MiSaTaQuWa formulation. Detweiler & Whitin g [3 provided an alternate 
formulation based on geodesies of a perturbed spacetime. These were followed by |15l - [l9| . which developed more 
practical methods for computing the actual self-force in Schwarzschild and Kerr spacetimes. They employed the 
so-called "mode-sum scheme" in which the scalar, vector or the tensor perturbation is decomposed in terms of 
corresponding spherical harmonics. In the case of GSF, a tensor spherical harmonic decomposition of the retarded 
metric perturbation ft,^,y(i,x) is performed. Then, the resulting 10 second order coupled partial differential equations 
are solved numerically at each tensor mode {£, m). The resulting metric fields and their derivatives are added together 
in certain combinations. These combinations are then translated from tensor modes to scalar {I, m) modes to yield 
individual I modes of the 'full' GSF given by Eg. ([7^. As the full GSF is singular at the location of the particle, 
a regularization procedure is undertaken. In MiSaTaQuWa formulation, this is done by decomposing the divergent 
'direct' part of the GSF into scalar spherical harmonics then removing these from the full GSF at each / mode. The 
resulting regularized / modes are finite and yield a convergent sum. This sidesteps the issue of dealing with infinities. 
The final GSF is then given by summing over the I modes from zero to infinity. 

The mode-sum scheme has thus far been implemented by several groups for SF computations [2l| - (3l| . Most of 
these have been for scalar field SF or looked at simplified cases for GSF computations (in Schwarzschild) such as 
radial infall or a static particle. The GSF for circular orbits in Schwarzschild was first successfully calculated (in 
time domain) by Barack & Sago [32| . This was soon-after followed by independent calculations by Detweiler [s^] and 
Berndtson [33|- Although these calculations used different gauges and methods, by comparing the effects of the GSF 
on gauge invariant quantities derived by Detweiler (s^j, these three independent GSF computations were shown to 
be equivalent (sH, |3g|- The state of the art for GSF computations is the recent work of Barack & Sago on eccentric 
orbits in Schwarzschild spacetime [l|. Some progress has also been made for GSF computations in Kerr spacetime, 
the state of the art being the work of Warburton & Barack [s^ on scalar field SF for bound (eccentric, equatorial) 
orbits in Kerr spacetime. This work was successfully implemented in frequency domain using the recently developed 
method of extended homogeneous solutions [s^. This was a very important step in the efforts to compute the GSF 
using frequency-domain methods. The method of extended homogeneous solutions successfully avoids the 'Gibbs 
phenomenon' that causes the radial derivatives of the metric fields to be averaged out across the point particle as 
opposed to displaying the expected finite jump there, which is the result of modeling the particle as a delta-function 
distribution. A very thorough introduction to the fundamentals of the self- force problem is presented by Poisson . 
In addition, a recent article by Barack [4l| overviews the current state of the field. 

Our aim in this part I of the series is to provide a fast framework for computing the GSF that can be used for orbital 
evolutions. For this reason, we have chosen to work in frequency domain (f-domain). Berndtson [35| was the first to 
successfully compute the GSF for circular orbits in Schwarzschild using f-domain methods, but his method differs from 
ours and his work is unpublished. Starting with Regge & Wheeler's (RW) standard tensor harmonic decomposition of 
the metric perturbation |42| . Berndtson solved the field equations in Lorenz gauge by relating the gauge invariant RW, 
Zerilli master functions [43| to the unknown metric fields of Lorenz gauge. It turned out, however, that he did not 
have the correct expression for the contribution of the monopolc mode to the GSF. But when he adopted Detweiler 
& Poisson's [1^ solution for this mode, the results he obtained for the GSF matched those of [s^. His results also 
highlighted the key advantages of a f-domain computation, namely, higher accuracy and faster runtimes compared to 
time-domain methods. 

Despite the evident success of Berndtson's approach, it is our feeling that our f-domain approach is better suited 
for extension to Kerr in that it relies less on the spherical symmetry of the background spacetime. As there currently 
exist no tensor spheroidal harmonics, we must rely on a tensor spherical harmonic decomposition of the metric 
perturbation in Kerr. The problem then is that the resulting ordinary differential equations (ODEs) couple between 
different multiple modes, not just metric fields. However, the principal parts of the ODEs remain uncoupled and it 
is possible to numerically solve the resulting system of coupled ODEs by treating the extra couplings as new source 
terms. We refrain from elaborating further as this problem is beyond the scope of this article but our longterm 
research program includes tackling these issues. 

The obvious advantage of working in the f-domain is that one deals only with ODEs, which can be solved efficiently 
using numerical methods. Furthermore, in f-domain, there are no instabilities associated with the non-radiative modes 
(monopolc, dipolc) that one encounters in the time domain [30l [s^. However, there are downsides to working in the 
f-domain. One is that f-domain methods work only for bound orbits. Also, it is generally thought that f-domain 
computations of GSF are intractable beyond eccentricities of approximately 0.7 [4J|. The breakdown of f-domain 
computations is caused by the fact that as the eccentricity increases, there are more and more radial frequency modes 
per given azimuthal mode. This significantly augments the runtimes of numerical computations. Eventually, one 
expects to reach a threshold eccentricity at which the use of time-domain methods becomes numerically more efficient. 
It is likely that f-domain methods become computationally inefficient (compared to time-domain) at eccentricities 
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higher than 0.7. We hope to empiricahy determine this threshold value in our future work. However, this may not 
necessarily present a problem since EMRI orbits circularize [i^ as they shrink toward the last stable orbit and despite 
recent findings that we should expect to see EMRIs with moderate eccentricities in the LISA bandwidth, 

most of the eccentricity will have been reduced by the time the compact object begins its final year of inspiral so that 
there should be plenty of EMRIs with eccentricities < 0.2 for LISA to detect. For such eccentricities, we expect an 
f-domain code to be significantly faster than its time-domain counterparts. 

As the GSF is a gauge dependent quantity (as is the orbital radius), we must address the issue of gauge choice used 
in our GSF computations. Lorcnz gauge is a common choice in pcrturbative studies of curved spacetimcs at linear 
order. One is motivated by this gauge choice because it retains the local isotropy of the delta-function singularity 
used to model the compact object j4a |. It also casts the field equations in a fully hyperbolic form, which is suitable 
for time-domain calculations. On the other hand, the perturbed field equations are generally more tractable in gauges 
like the Regge- Wheeler (RW) [i^ or the radiation gauges 47 1. However, thanks to the work of Barack & Lousto [2G |. 
we now have access to all of the field equations in Lorenz gauge and can follow an "all-Lorenz-gauge" path. This 
is especially desirable in the mode-sum scheme because the multipole modes of the metric perturbation {hj]J){t,r) ) 
in Lorenz gauge are continuous (C°) at the location of the particle. This is not the case, for example, in RW gauge 
where the source contains a derivative-of-delta-function term in addition to the usual delta function. Therefore, the 
so-called "master functions" used in the RW formalism exhibit a jump-discontinuity (C~^) at the location of the 
compact object. Finally, one can compute only the radiative {£ > 2) modes of the perturbation using approaches 
based on RW gauge 



Our treatment here is mostly based on the work of Barack & Lousto [20| (henceforth BL) and Barack & Sago 
2007 dH (BS), which use the mode-sum method in Lorenz gauge. We begin with the linearized Einstein equations 
in Schwarzschild background in Lorenz gauge. We then rewrite the field equations using tensor spherical harmonic 
decomposition of the metric perturbations. This decouples the angular part of the field equations. The resulting 
set of 10 second order partial differential equations are separated into 7 even and 3 odd parity equations. Next, 
we go into the frequency domain and obtain 7 © 3 second order ODEs. For a generic bound orbit, we would need 
to sum over radial and azimuthal frequency modes to work in f-domain, but for circular orbits we have only one 
fundamental (azimuthal) frequency. Therefore, the crucial step in moving to an f-domain computation for circular 
orbits is supplying appropriate boundary conditions for the metric fields. Here, we present these boundary conditions 
(BC) for the first time. 

With the BC specified, we numerically solve the coupled homogeneous ODEs then impose junction conditions at 
the location of the particle to construct the inhomogeneous solutions. Once we construct all the metric perturbations 
and their derivatives at the particle, we compute the GSF by using the formulae derived in BS. This gives us what 
is called the "full self-force" . It contains a 'tail' contribution, which we interpret as the relevant physical piece and 
a 'direct' part, which must be removed via the appropriate regularization procedure. It should be iterated that the 
initial decomposition of the metric perturbation is done in tensor spherical harmonics, whereas the regularization is 
performed using scalar spherical harmonics. This requires us to translate each tensor (£, m) mode to various scalar 
{l,m) modes before regularizing. This causes a single scalar mode / to couple to many tensor modes £. The formulae 
for these couplings have been derived by BS. Here, we use their results to compute the GSF. 

For circular orbits, only the r-component of the GSF needs be regularized. In the mode-sum scheme, this is done 
modc-by-mode at each scalar multipole I where the singular piece is decomposed in scalar spherical harmonics then is 
removed from the full self-force at each I. The resulting regularized I modes have large-^ behavior, which yields a 
convergent (albeit somewhat slow) sum over I. The physical self-force is obtained by summing over all the individual 
regularized I modes and finally adding a large-^ "tail" that estimates the total contribution due to I > Imax modes 
where Imax is the largest mode at which we actually compute the metric perturbations. 

Section [H] presents the field equations and their decomposition under tensor spherical harmonics. In section Hill we 
go into f-domain by Fourier transforming the time dependence of the metric fields in azimuthal frequency modes. We 
then separate the resulting field equations under their parity and calculate the BC for each case separately. Once the 
BC are known, the numerical ODE solver integrates the field equations to yield the homogeneous solutions. Using 
these, we assemble the inhomogeneous solutions, which we use in section ITVl to construct the full GSF, which we then 
regularize. Finally, we compute the tail contribution to the r-component of the GSF. The results are all displayed 
in section IVl where we compare the t-,r-components of the GSF computed by our code with that of BS. We find an 
excellent agreement with BS within their error bars for orbital radii up to ~ lOOGM/c^. 

Throughout this article, we use geometrized units with G = c = 1. = {t,r, 9, 0) are the standard Schwarzschild 
coordinates and r denotes proper time. We follow the usual convention of (— , -I-) for the metric signature. 
Finally, owing to the spherical symmetry of Schwarzschild spacetime, we work with equatorial 6 = 7r/2 orbits without 
loss of generality. 
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II. FIELD EQUATIONS 



The physical set-up is that of a point particle with mass /i in a circular orbit with radius tq around a Schwarzschild 
black hole with mass M. The particle interacts with its own gravitational field and thus feels a net force which moves 
it off the geodesies of the background spacetime. The equation of motion for the particle in this context is given by 

fiu^W^u'' = Fl^sF, (1) 

where u^^ = dx^/dr denotes the 4-velocity of the particle, r is proper time, is the covariant gradient operator 
associated with the background Schwarzschild metric and i^Qgp is the gravitational SF. Imposing the condition that 
the 4-velocity remain normalized along the worldline i.e. m^m^ = — 1 on Eq. ([T|), we get the orthogonality condition 
on the self- force: u^F^^p = 0. For circular orbits, -^gsf'-^gsf calculated independently using energy balance 

arguments (49j because they are purely dissipative. However, in the case of eccentric orbits all non-zero components 
of the SF will be made up of both dissipative and conservative parts. The orthogonality condition is useful because 
it gives us a simple way to obtain one out of the three components of the GSF (fourth component Fgp = because 
of spherical symmetry). 

To obtain the GSF in this "forced geodesic" picture, we must solve the perturbed Einstein's equation in a non-flat 
background. Schematically, the field equations have the following form 

GtiA9tii' + hfii,] = SttT^^, (2) 

where G^i, is the Einstein tensor, which is a functional of the spacetime metric g^^ ~ g^^ + h^^i, and T^^ is the 
energy-momentum tensor sourced by the point particle. Here, g^jy denotes the background (vacuum) Schwarzschild 
metric and /i^i/ is the perturbation due to the point particle. As is standard with current GSF computations, we retain 
only the linear order ©(u) perturbation. There are ongoing efforts to incorporate second order perturbations in the 
calculations of GSF [50l.[5l|. but the current formulations are not yet ready for use in mode-sum GSF computations. 

After keeping up to 0{h^^) terms in Eq.®, we substitute G[g\ = into Eq.® since g^^ is the metric of a vacuum 
spacetime. Wc make two more simplifications, which are standard: first, we change from using /i^j^ to the trace- 
reversed hfj,i, via h^u ~ h^i, — \g^uh. Then, we pick a gauge. For reasons explained above and detailed in the cited 
articles, we choose to work in Lorenz gauge where S/^h'^" — 0. With these modifications inserted into Eq. ([2]) we 
obtain 

nh^, + 2irfjKp = -i6^T^., (3) 

where □ — V^V^. The energy-momentum tensor is given by 

T^u = fi {-g)-^/^S^[x''-x^{T)]u^u,dT, (4) 



where x^{t) denotes the position of the particle. The proper time r is related to the coordinate time t via dr ~ 
(v})~^dt. Finally, g is the determinant of the Schwarzschild metric equaling — Tq for 6 ~ Tr/2. 

As it stands, Eq. ([3]) represents 10 coupled 2""^ order, partial differential equations (PDEs). We can simplify these 
by separating out the angular part. To this end, we decompose hfi^(t,r) using tensor spherical harmonics, which 
form a 10-dimensional basis for any rank two, symmetric 4-dimensional tensor field. The components of the metric 
perturbation are decomposed as follows: 

10 

V(t,r) = ^5]5]/iW^™(t,r)r«^™(0,0;r). (5) 

i,m 2—1 

The explicit expressions for Y^J are presented in BL. We modify them slightly here: Y^J^^^^ = oS^' Y^J where 
a*-*'^ constant coefficients defined in BL. Now angular variables decouple and the field equations become (at each m) 

+ M g5 /i^^'^'" = . (6) 

where / = /(r) = 1 — and n^c is the usual scalar field wave operator: 



sc ^ 



(7) 
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are the source terms obtained from decomposing T^i^ in tensor spherical harmonics. They are given by 



where Eq ^ {1 — 2M /ro)/y^l — 3 A// ro is the dimensionless energy of a test particle (/i = 0) on a circular geodesic 
with radius j-q. Given the orbital angular frequency f^o = d4>Q/dt ~ (M/vq)^^^ , the constants a^*-* are: 

= /oVro, a'^' = /o/ro, a^"' = ronl[£{£ + 1) - 2m% 
a(2).(5).(9) = 0, =2/ol7o, 

a^'^ = ro^l (9) 

where /o = 1 — 2M/ro. Note that the (i) = 2,5,9 equations are sourceless. The spherical harmonics are given by the 
usual formula 



r^™(6», (f>) = c^™P^™(6l)e™'*. (10) 

ilyno 

Eq. ([5]) using the following expression: 



P^™ (9) are the associated Legendre polynomials and c^™ = y 2£±i . "We can further rewrite the second line in 



dg[Y'^{7:/2, 0o)]* = [^C,+i^„c,+i,™ F^+i''"(7r/2) - {£ + 1)Cm^ P^-1'™(^/2)] e"™^" 

= Jft'^'^'e-™'^" (11) 



where Cfo„ - y (2f+i)(2f-i) • 

The l^.^^ft,'^-') in Eq. (j6]) contain the coupling terms between different field equations. In the next section, we will 
show that up to 5 field equations couple together for certain modes, but things will not get any more entwined than 
that. The expressions for AA^'^^-^hS^^ are lengthy and have been given in detail in [l[, [2^ and [s^l so we omit them 
here. We will however present the field equations in frequency domain in section ITlll 

Eq. (|5]) substituted in to Eq. ([5]) gives us the Einstein field equations in their simplest form that we can reach in 
Lorenz gauge. From this point, one can either go into time domain and tackle the problem of solving these coupled 
PDEs or one can go into frequency domain and deal with ODEs that require boundary conditions. In the next section, 
we solve the field equations in frequency domain in Lorenz gauge for the first time. 



III. FREQUENCY-DOMAIN SOLUTIONS OF THE FIELD EQUATIONS 

Here, we begin by decomposing the metric fields /i*-*' (t, r) into frequency modes. In the case of circular orbits, there 
is only one frequency: i7o = (M/rg)^/^. So the harmonics of circular motion are given by = mflo- For elliptical 
orbits, the frequency modes will be a combination of azimuthal and radial fundamental frequencies: LJmn = mfl^+nQj.- 
For circular orbits, metric fields are decomposed as follows: 

h^^)'"^{t,r) = where = mQo- (12) 

This reduces the 2-dimensional hyperbolic equations ((6)) to a set of 2"^^ order, coupled ODEs, which can be numerically 
solved much more quickly than PDEs encountered in time-domain approaches. In the case of a scalar field in 
Schwarzschild spacetime, the problem in f-domain reduces to a single inhomogeneous ODE. The standard procedure is 
to numerically solve for the homogeneous inner (r < rg) and outer (r > tq) solutions then construct the inhomogeneous 
solution by imposing the correct junction conditions at r = rg. For the computation of the GSF, the same procedure 
applies but now for many coupled fields, some of which have delta-function sources and others no sources at all. In 
section llll A 11 we explicitly show how we construct the inhomogeneous solutions from coupled homogeneous solutions. 
The system of 10 coupled, second order homogeneous ODEs can be written as 

^^^^ - mm{r)4lir) 4.M« i?« (r) = 0. (13) 
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where is the Reggc- Wheeler tortoise coordinate with dr^/dr = f ^, is the Fourier transformed version of 



A^(^), and 



Vim{r) 



2Mf , £{£+!)/ 



(14) 



The field equations (jl3|) are not all coupled to each other; our 10-dimensional basis splits under parity very much like in 
Regge- Wheeler gauge. The (i) = 1, . . . , 7 basis elements of the tensor spherical harmonics are even and the {i) = 8, 9, 10 
basis elements are odd under parity transformations. For circular orbits, even, odd mean that £ + m = even, odd. 
Eqs. now decouple completely under these two parity sectors so they can be solved completely independently. 
Furthermore, because the spherical harmonics in the source terms ^ give [Y^"^{tt/2, (j)o)]* = for £ + m = odd and 
d0[Y^™{TT/2, (po)]* = for £ + m = even, the odd parity solutions are trivially zero for an even mode and vice versa 
for even parity solutions. That is rW---('<"> for £ + m = odd and = for £ + to = even. 

Similarly, the four gauge equations coming from the Lorenz gauge condition V^/i'^j^ = also decouple under parity 
with three equations falling under the even parity sector, leaving only one for the odd sector. The gauge equations at 
each {£, m)-mode are 



(1) 



/ 



(2) 



0, 



(15) 



- - ^ (ri?(5) + 2i?(^) + £{£ + l)i?(6) - = 0, 



(16) 
(17) 



_ _ L (^i?(9) + 2i?(9) - = 0. (18) 

Here and henceforth, we omit writing the modal indices £, m as well as the functional dependence on (or r) for the 
sake of brevity. It should be assumed that each field equation presented holds for a given £, m mode unless stated 
otherwise. 

Thanks to the gauge equations, it turns out that not all the even (or odd) equations need to be solved simultaneously. 
As we have four gauge conditions, we have only 10 — 4 = 6 degrees of freedom (d.o.f). These split as 4 + 2 under 
parity. But because of the particular form of the field equations in the even sector, we must solve 5 coupled ODEs 
together, construct the inhomogeneous solutions then use two gauge equations to obtain the fields i?^^^ and i?*-^^ (more 
on this later in section UlIBp . In the odd sector, we solve the two coupled (i) = 9, 10 equations together then use 
the odd gauge equation to obtain R'-^K This procedure of solving the equations in stages is called "the hierarchical 
solving scheme" by BL. It involves first numerically solving only the ODEs that couple to each other then using gauge 
equations psp - (jlSp to determine the remaining unknown radial fields. The number of equations one has to solve 
changes depending on the values of £ and m. For a generic even mode {£ > 2, to > 1), one solves 5 coupled ODEs 
then uses two gauge equations whereas for a generic odd mode {£ > 2,to > 1), only two coupled ODEs are solved 
numerically then one gauge equation is used. There are also non-generic modes such as the monopole (^ = 0); the 
even, odd dipoles (^ = 1,to = 1,0) and the static (to = 0) even, odd modes. Analytic solutions have been explicitly 
provided in [33| for the monopole, and by BL for the odd static modes. The even dipole {£ = l,m = 1) and the 
£ — even static modes are solved numerically, but have fewer number of non-zero fields. We present all the different 
cases for both even and odd parity sectors and the hierarchical scheme for solving the field equations in table |T] 



A. Odd Sector 



We begin with what we call generic odd modes (to > 0). We will consider the static odd modes (m = 0) later 
in a special subsection. As explained in the hierarchical scheme, here we solve the coupled (i) ~ 9, 10 equations 



7 





Even {£ + m = 2N) 


Odd H 


' + m = 2Ar+ 1) 


e = o 


(i) = l,3,6^2 (A) 




no field 


i = i 


m = 1 : (i) = 1,3,5,6 2,4 


m = 


(i) = 8 only (A) 


e>2 


(i) = 1,3,5,6,7 ^ 2,4 




= 9, 10->8 




m = : (i) = 1,3,5 ^ 6,7 


m = 


(i) = 8 only (A) 



TABLE I: The hierarchical solving scheme for the ten field equations. The arrows — >■ indicate that we use the gauge equations 
to obtain the field to the right of the arrow. (A) indicates that the solutions are obtained analytically and N £'N 



together to determine i?'^' and R^-^^^ then use these solutions in the odd gauge equation to solve for i?'^^ The 
two homogeneous, odd parity field equations are 

where \ — {£ + 2){£ — 1). In order to get the correct numerical solutions to Eqs. ([T^ and (|20p . we must specify 
appropriate boundary conditions for the numerical ODE integrator. The boundaries are located on the event horizon 
(r — 2M) and at radial infinity (r = oo), which translate to r* = ~oo and r, — oo, respectively. A quick inspection of 
the structure of the ODEs (i) = 9, 10 reveals that as r, — ?> oo and r — > 2A'/(r, — >■ — oo), the term dominates in 
the potential and the ODEs ([TO)) and (PU)) asymptotically turn into standard wave equations. Thus, for the solutions 
at infinity and on the event horizon, we have the usual outgoing and ingoing wave behavior, respectively. Denoting the 
outgoing/ingoing homogeneous solutions by R'^ and R~ , respectively, we write the following ansatz for the boundary 
conditions: 

oo fc 

e-"'-* J2 (21) 

oo 

^-'^"-'■*E^9'io(^-2M)'=. (22) 

fe=0 

Clearly at r = 2M and r = oo we get the proper wave- like behavior. We must also specify dRf /dr [i = 9, 10) at the 
boundary points. Our numerical code uses as the integration variable so we actually need dR^ /dr^, = fdR^/dr 
for the BC. 

Numerically, we can not use infinities for the boundary points. For our code, we pick a range of e [— 65il/, — 55M] 
for the inner boundary, = — 65A/. which corresponds to r/M « (2 + 10~^^) is about as far 'in' as we can go due 
to double floating point machine accuracy. The choice for the outer boundary point Tout depends on £ and ujm as 
we demand that the outer boundary be located in the wave zone, which translates to Tout 3> {£ro)/u>m- So we opt 
for an adaptive outer boundary at each {£, m) where rout = 50 (£ro)/wm- The ratio of 50 was chosen after numerical 
experimentation. Larger ratios mean larger runtimes for the computation of the homogeneous fields, and smaller 
ratios call for more terms in the series in Eqs. (PT|) . (|22p for numerical convergence. 

Note that the sums for the BC in Egs. pTj) . ([22)1 are infinite. However, because we solve the coupled field equations 
numerically, we must truncate the sums at some k = fcmax- We numerically determine this fcmax for each of the sums at 
every {£,m) such that the next term in the summation has absolute magnitude less than 10~^^. We also numerically 
check that each sum converges. 

The coefficients and b\. are unknown and must be determined by substituting our ansatz into the field equations 
then constructing recursion relations for the fc"^ coefficients a], and b]. out of a^,^j, and h\,^i.. The recursion relations 
for the outer BC for R^ and i?^g are as follows: 

2iujk 4 ^ Ck-i al_^ + Dk-2 al_.2 + Eks a^a + 2^-1 - WMal% + 12M^al%, 

(23) 

2zLok 40 = af_i + Jk-2 aj}l2 + Kk-3 al'L, + 2Xal„, - UlXal.^, (24) 



/ / 45M 
V r 



/ ^^(lo)_^^(9)^ 



2r2 



(19) 
(20) 



R. 



•9,10 
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where 

Ck = U'liojk + fc(fc + 1) - L - 4, 4 = AMiiok + k{k + l) - L + 2, 

Dk = -6Mk - iMk^ + 24A/ + 2ML, Jk = ~6Mk - AMk^ - 6A/ + 2ML, 

Ek = 4A/2(fc^ + 2fc-8), Kk^'iM'^ik^ + 2k + l). 

Here and in all other recursion relations that we present, lo denotes Um = mU,^ and L = i(£ + 1). The recursion 
relations are rather cumbersome, which is why we will refrain from presenting the rest of them in the main body of the 
paper unless we refer to them directly (as done in section llIIB4p . All the recursion relations are listed in appendix 

El 

The recursion relations must be started off by specifying the values for the leading terms. In the case of odd parity 
equations, these first terms are Oq'^'^ and ^q'^" with o^'^q = and b^i.<S) ^ ^' This gives us 4 free parameters to specify 
every time we wish to solve the system of coupled ODEs. Since we have one gauge equation and three field equations 
for both inner and outer homogeneous solutions, we end up with 2 x (3 — 1) = 4 degrees of freedom. These d.o.f. 
are manifest in our freedom for choosing the values for Cq'^'^ and &q'^°. In the next subsection, we show how to pick 
suitable values for these coefficients and construct the inhomogcncous solutions. 

The final remark concerns the nature of the BC specified above. As can be clearly seen, the ingoing/outgoing wave 
conditions for the BC yield complex numbers. Therefore, we must construct complex solutions for the homogeneous 
fields Rf . A quick inspection reveals that the real and imaginary part of the complex fields i?*^'^ completely decouple 
in the field equations ((T9| and pO)). As a result, we simply solve each given ODE twice: once with the real part of the 
BC and a 2"=^ time using the imaginary part of the BC. Wc then combine the two numerical homogeneous solutions 
under one complex solution that wc also call Rf . Recall that wc already have to solve the homogeneous ODEs twice 
to get the inner (— ) and outer (+) solutions and now twice more for the real and imaginary parts. In total, at each 
generic odd mode, wc must numerically solve the system of coupled ODEs 2x4 = 8 times. 



1. Obtaining The Inhomogeneous Solutions 

To obtain the true, inhomogeneous solutions — which are sourced by 5-functions — we must impose junction 
conditions on the coupled homogeneous solutions. Recalling that the inhomogeneous solutions i?^*) must be C° fields, 
the two conditions are continuity at tq and the correct jump of dR'^^^ /dr across rp. Because we have coupled fields, 
we must construct the inhomogeneous solutions from linear combinations of homogeneous solutions. We use standard 
methods of constructing a linearly independent basis of homogeneous solutions and imposing the correct junction 
conditions to assemble the inhomogeneous fields. Below, we briefiy outline this procedure. 

As mentioned before, in the odd sector we have a total of 4 d.o.f. so we construct a 4-dimensional basis from the 
homogeneous solutions Rf and Rf^. We do this by exploiting the freedom we have in choosing the initial values 
for the coefficients a)^'^g , 6^.'^g that start the recursion relations (|23p - (p4)) . A linearly independent 4-dimensional 
basis can be constructed for the homogeneous solutions i?9.io by setting {a^^a}^) = (1,0) then (0, 1) and the same 
for (6g,6j'^). These determine our basis vectors at the point of interest, namely r — tq. We label these solutions by 
i?g^'^, and R!i^^^ , R!\f^ . For example, i?g^'^,_R^g^ are obtained by setting ag = 1 and a}^ ~ then solving the 

coupled ODEs ([19]) and for Rf{rQ) and i?]*g(ro). Recall that since the boundary conditions are complex, the 
basis vectors are complex as well. Finally, we follow the same procedure for the r-derivatives. We label the inner and 
outer basis elements for the r-derivatives Sr^g^'^, <9ri?io^ and drR^^g^^ ^drR^^il"^ ■ i'^ ^^is notation, drP^il~ stands 
for di?]~o/dr|ro obtained by setting 6o = and h}^ ^1. 

We label the inhomogeneous solutions by , i^^^^ . These are constructed from ^9^10,^9^10 respectively. The 
inhomogeneous solutions are obtained by imposing the standard junction conditions: (1) Continuity at tq: i?|*^(ro) — 
^out(^o), (2) The following jump for the r-derivatives at ro: 



dRolt 



dr 



dRf^ 



dr 



i^^^f^x^°-.J«, (z) = 9,10, (25) 
/o 



where J°'^'^ is given by Eq. ([TT|) . To impose these conditions for our basis of homogeneous solutions, we form a 
4x4 complex matrix containing the fields i?!"''^, Sri^P'^ listed above. The inhomogeneous solutions R^^o^t^ are 
constructed from linear combinations of the homogenous solutions multiplied by unknown complex coefficients Xj . To 
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determine these coefficients, we must solve tlie following matrix equation: 



[1]- 



— drRg 



r[2]- 
— -Kg 

n[2]- 
-^10 

— drRg 

— drR!?r 



tig 

drRg 



n[2] + 

tig 

n.[2]-' 
^10 



drR 



[1]+ 

10 



drR. 
drR 



2] + 
3 

2] + 
10 



\ 






( ' \ 





























(26) 



The right hand side (RHS) of Eq.(p6| ensures the continuity of the inhomogeneous solutions and imposes the correct 
jump value J^*^ on the first derivatives. Recall that because a^^^ = (see Eq.®), we have j'^^ — 0. We solve for 
the complex xi, . . . , X4 by using standard numerical matrix inversion algorithms. Once we know the cci, . . . , X/^, we 
construct the inhomogeneous solutions at the location of the particle. These are given by 



dR^^ 



x,R^- + x^R^- ^ x,R^' 



X4R, 



[2] 



dr 



= XidrRj 



X2drR. 



[1]- 



dR. 



dr 



X3drR^^^^ + XidrR^i^^^ , 



(27) 
(28) 



where i = 9, 10. Although the continuity of i?^^-', and dR'^^^ /dr (because J^^-* = 0) is analytically exact, because 
the coupled ODEs are solved numerically, we will inevitably have a small violation of continuity at r = rp. This 
is caused by the numerical matrix inversion. Usually, the numerical inversion algorithms are very robust and the 
discontinuity in the fields is ~ 10~^^ — 10~^* for most modes. However, for a few special modes, this error becomes 
much more significant. We will comment more on this issue later in section flVBI 

We take the solutions pT]) and substitute them into the odd gauge equation p8)) to solve for i?(®)(ro). After this 
step, we obtain dR^^^ /dr a.t r — tq by differentiating the gauge equation ((T5)) with respect to r and using the field 
equation to substitute for d^R'-^^ term in /dr. RecaU that i?^^) j^^s a non-zero (5-function source thus it 
exhibits the standard jump discontinuity at tq given by Eq. (j25p . Therefore, we must compute dR^^^ /dr\rg twice: once 
as r — >■ Tfj" then again for r ^ r^ . Since i?'^-'(ro) and its ± r-derivatives are obtained algebraically from Eq. lfTSl) — by 
inserting the numerical solutions i?*^^^''^*^^ (rg), dR^^''''-^^'' /dr\rg — we expect the error in the continuity of i?'^-'(ro) to 
be comparable to errors found for i?(9)'(i°)(ro). Indeed, we find that the offset in the continuity of R'-^^rg) is ^ 10 
Similarly, the relative error between J*^^) and the jump of dR^^'^ / dr\rg is ^ 10~^^. 

As mentioned above, we have to solve the set of coupled ODEs 8 times for each odd parity mode: twice owing to 
the fact the BC are complex, and 4 times because we construct the inhomogeneous solutions from a 4-dimensional 
basis of homogeneous solutions. Doing a run up to e.g. imax = 18, we end up with 81 generic odd modes, which yield 
a total of 81 x 8 = 648 times that the coupled set of odd ODEs must be solved numerically. 



2. The Static (m = 0) Odd Modes 
As shown in BL, the rn ~ odd modes have analytic solutions. Since 

J(10) 

cc m ~ and J(9) = 0, we trivially have 
that i?'^-* = = for these modes. Therefore, we solve a sing le ODE for R.'-^l For the case of £ = 1, the ODE 

simplifies to a well known form, which has the following analytic solution: 



1 o j (r/rof, r <ro 

>3 (ro/r), r > rg. 



(29) 



where l3e=i = leVs^fQ^Eo^o . For £ > 1, the inner (r < rg) homogeneous solutions exhibit the standard power law 
behavior: ~ r^+^. As for the outer solutions (r > rg), we have something that is of the form r~^(l + In/). These 
scale as r^^ as r — > od, which is regular. The details of how these analytic solutions are constructed are given in 
section IIIC of BL, which is why we refrain from elaborating more here. We also omit the explicit expressions for 
these static, £ > 1 solutions in this article. The interested reader should peruse BL ([l^l)- In summary, the overall 
static, odd solutions are given by — restoring the modal indices — hi^y^ = /i(io)^o = Q and the non-zero fields 
which are constructed analytically . 



B. Even Sector 



For the generic, non-static case of even modes, we have 7 field and 3 gauge equations thus a total of 2 x (7 — 3) = 
d.o.f. However, an inspection of the even parity field equations as they are written in Lorenz gauge ([l|, [20j . (32 
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reveals that we must simultaneously solve 5, not 4, coupled ODEs. As before, we numerically solve the homogeneous 
ODEs then obtain t he inho mogeneous solutions by employing the standard techniques for coupled fields, which we 
illustrated in section UlI A II The 5 homogeneous coupled ODEs in the even sector are the (i) = 1, 3, 5, 6, 7 equations 
written in the following form; 



r J 



R 



(5) 



fR 



(3) 



^ ri - ^) R^' 



dlR<-'^ 



_ ^_ 



m^R 



4.5Af 



(30) 
(31) 

(32) 
(33) 

(34) 



In this article, we follow the convention of BS ^3 for the field h'^^\ which is different from that of BL VK 



fiil/f. Recall that / = 1 - 2M/r, L = £{£ + 1), A = (^ + 2){e - 1) and Vi„r is given by Eq. ([H 



here 

Next, we must specify the boundary conditions. As was the case with the odd sector fields, we impose the same 
ingoing/outgoing wave conditions on the event horizon and at radial infinity, respectively. We once again use R~ , Rf 
to denote the ingoing, outgoing homogeneous solutions, respectively. For the inner/outer BC, we use the same ansatz 
as before 



R7 = e 



R 



OO A 



(35) 
(36) 

34l) to derive new recursion 



for i = 1,3,5,6,7. Once again, we substitute these ansatz into the field equations ((30 

relations for the coefficients a^, 6^ in Eqs. ([55]) . ([55]) . The sums are of course infinite but we truncate them at some 
k — fcniax as we did before. The recursion relations for the outer coefficients and inner coefficients h\. are given in 
appendix Rl 

With the coefficients , \)\, determined, there still remains one critical issue that pertains to the total number of 
degrees of freedom to use: in the even sector, we have 5 ODEs that can not be decoupled from each other, so we 
must solve all five simultaneously, but we have 8 d.o.f in the even sector, not 2 x 5 = 10. So, there must be an extra 
condition on each set of 5 BC for inner and outer homogeneous solutions. For the outer solutions, this extra condition 
is a constraint on the coefficients a^, which is given by the even gauge equations: 

0^ = 0. 



(37) 



We repeat this procedure of eliminating the 5**^ degree of freedom from the inner homogeneous solutions by making 
use of the gauge equations. After some manipulation, we reach the following condition on the coefficients 63: 



[( + 1) + 4Ma;(l - 4Mia; + £(^ + 1) ) ) 6j + i (1 + IGM^a;^)^^] 
2Muj{\ + 16M2w2) 



(38) 



So all of the coefficients h\ are entirely determined from h\, h% and the recursion relation (jA22p . With the conditions 
([57)1 and imposed, we are left with the expected 8 d.o.f. 

Eqs. (|37|) and ((38)) tell us that our 8-dimensional basis of inner and outer homogeneous solutions is constructed by 
using the recursion relations (jA6|) - (|A29|) for the BC with 69, ^q, 69} and {aj, Og, Oq, Og} as the sets containing 
the 8 free parameters for the inner and outer homogeneous solutions. We construct our basis of linearly independent 
homogeneous solutions by numerically determining the basis vectors that span the solution space. Each basis vector of 
the outer homogeneous solution space is obtained by setting one of the coefficients {aj, Oq, ag, Oq} equal to 1 while the 
other 3 equal 0. We do this a total of 4 times, e.g. {aj, ag, a%, a^} = {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0} and {0, 0, 0, 1}. 
This procedure is repeated with 6q, 6q, 5q} for the inner solutions. This yields 8 basis vectors for constructing 
the 8-dimensional linearly independent homogeneous solution space. Given that the system of ODEs must be solved 
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twice because the BC are complex, we reach a total of 8 x 2 = 16 for the number of times we must numerically solve 
the field equations at each even mode. For example, for £ running up to 18, we have a total of 89 generic even modes, 
which means that the coupled ODEs arc numerically integrated a total of 89 x 16 ~ 1424 times. This is what takes 
up the main bulk of our numerical computation time. We will say more about this later. Next, we construct the 
inhomogeneous solutions. 



Inhomogeneous Solutions 



In subsection IIII A 11 we showed in detail how to construct the inhomogeneous solutions from the inner and outer 
homogeneous solutions. Here we do the same with the even parity solutions. Our basis of homogeneous solutions 
is now 8-dimensional and is spanned by K^,drR^ with i = 1,5,6,7. In accordance with the notation of subsection 
IIII A 1[ we label the basis vectors (the homogeneous fields i?f 5.6 7) by P^f^- Similarly, for the derivatives, we use 

drP^^^ ■ For example, iJ^^'"*" stands for obtained by setting a\ = \ and Oq = Oq = Oq = and S^iig"^'" is dR^ /dr 
with bl = l and 6j = 5^ = 6^ = 0. 

To construct the inhomogeneous solutions, we impose the junction conditions on the homogeneous fields and their 
r-derivatives in the form of an 8-dimensional complex matrix equation: 





(39) 



O4XI is a 4 X 1 array of zeros imposing the condition of continuity for the inhomogeneous fields and 



J 



(40) 



The complex, inhomogeneous fields i?*-') at r = are given by 



(41) 



Similarly, for the r-derivatives of these fields at r = tq, we have 

dR^ 



dr 



dR. 



(0 



dr 



4 



(42) 
(43) 



We still need to determine the inhomogeneous field R''^^ and its r-derivative at tq. Recall that in order to form the 
linearly independent basis of homogeneous solutions we had to solve a system of 5 (not 4) coupled ODEs together. 
However, the homogeneous solutions R^ and their first derivatives R'^ arc not part of our basis because they are 
constructed from linear combinations of the other basis elements as shown in Eqs. ([57]) & (|55| . With the basis of 
homogenous solutions at hand, R^'^\rf)) , dR^'^'^ / dr\rQ are simply given by 



^:iut(^o) 


4 

= E 


x,Rf-^j: 










dr 




4 

-E 




ra 


i=i 




"-''■out 

dr 




4 

-E 


Xj-^-4^rR]^^~^ ■ 









(44) 
(45) 
(46) 
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The remaininff two fields and are extracted from tfie even parity gauge equations p6)) . pT)) . Their r- 
derivatives are obtained by differentiating these gauge equations with respect to r and substituting the relevant parts 
of the fields equations (i) = 1, 3, 5 for the R'^^\ R'^^\ R'^^^ terms that arise from r-derivatives of Eqs. ([T6)) . ()17p . 

Although (ro) = i?o^t(ro) analytically, because we invert the complex matrix numerically, we are bound to have 
small discontinuities at tq as we did with the odd parity fields. We checked the relative error in the continuity of the 
fields at r = To and found that it is at most 0(10 ^'^) for tq < lOOA/ and £ — m = small. However, 

we find that for tq > lOOAf , as £ — m — >■ 15, the violation of the continuity of the field i?'^) grows up to 0(10"^) in 
relative size. For i — m > 30, this violation climbs up to O(10~^). Clearly, for large orbital radii and large £ — m, 
the numerical matrix inversion becomes less accurate. A quick check of condition numbers c for the matrices in Eq. 
([39]) shows that c > 10^^ for the problematic cases mentioned here. We explain the cause of this in section HVBI 
However, it is only the field R^^^ that exhibits the bad discontinuities; the fields i?(i)'(6),(7)^ which also come directly 
out of the matrix inversion, have continuity violations that are consistently at least three or more orders of magnitude 
smaller. As expected, larger inversion errors persist in the fields R^-^\R^^'> (and their r-derivatives) because these 
are constructed from gauge equations containing i?*^^) and its first and second r-derivatives. As far as we can tell 
this matrix inversion error, which we quantify by the numerical discontinuity of the fields i?(2),(4),(5) at r = ro is our 
largest source of error. We will say more on this inversion error in section flVBI 



2. The Even Dipole (£ ^ l,m ^ I) Mode 

The even parity dipole mode is non-radiative (£ < 2) thus represents a shift in the orbital angular momentum, 
which can be interpreted as a rotation of spacetime around its center of mass. For £ = l,m = 1, A = Oas well 
as aC^) = JC^) = 0. This gives /i('^)"(t,r) = 0, which results in 4 coupled ODEs. The (i) = 1,3,6 equations ((30l) . 
([5T|). ([55]) do not contain any i?^^^ terms as such they remain unchanged, as do the recursion relations for a\'^'^, 
displayed in appendix I A II However the (i) = 5 equation p2p does contain a Ai?'^' term, which is now zero so we end 
up with new recursion relations for the inner and outer boundary conditions for Rf. These are given by Eqs. (jA3ip . 
(|A32|) in appendix IA21 

With Rf = 0, we have 2 x (6 — 3) = 6 degrees of freedom for our basis of homogeneous solutions. The basis 
vectors are constructed from the homogeneous solutions obtained by using the BC generated from the sets {b^, 6q, 6q} 
and {ao,aQ,ao}, respectively. The ODE integrator solves the coupled system a total of 2 x 6 = 12 times. To obtain 
the inhomogeneous solutions, we construct a 6 x 6 complex matrix very similar to the one in Eg. pQp . but without 

the homogeneous fields R^^^ , drR^^^ . We solve the resulting matrix equation to obtain the values for the complex 
amplitudes xi , . . . , xg , which in turn, give us the values of the inhomogeneous solutions and their first r-derivatives 
at tq. The equations for the inhomogeneous fields i?(i)'(3),(5),(6)(-y,^-| arc identical to Eq. (jH]), (|44l) with xt = xs = 0. 
The fields i?(2)X4)(^g) 

arc once again obtained from the gauge equations and (fT7|) with i?'^^ = 0. 



3. The Monopole 1^0 Mode 

This conservative, non-radiative £ ^ contribution to the metric perturbations represents a shift in the mass of the 
small particle across r = rn. For this mode, the field equations simplify enough that analytic solutions have been found 
by Detweiler & Poisson [33| . The only non-zero fields are = /i(3) = i?(3)^ /j(6) = ^(6)^ contribute only 

to the diagonal (scalar) components of /i^^. In section III.D of BL, the solutions are displayed explicitly in terms of 
the components of As with the other modes, these are C° with the usual jump in the r-derivative across rg. We 
omit writing the explicit solutions here and refer the interested reader to section HI.D for the details. The extra 
important step we mention here is the rewriting of these analytic solutions — written as components of h^^^ in BL — 
in terms of h'-'\ Although this seems like a backward step, it is necessary in order to properly follow the algorithm 
for computing the GSF. We will elaborate more on this procedure later in section HVl 

The formulae needed to transform htt,hrr,hgg, h^^ to }S^\TS^\}S^^ are as follows f j20j): 

2V^^- V {htt + fhrr) , (47) 
2V^^r'j{hu- fhrr) , (48) 

4x/¥/i"V"^/i00 = A^/^fi-^r-\sm0y^h^^. (49) 

Note that the expression for h^^'' here looks different from the one given by BL in [2^. The reader may recall that this 
is because we use the ft-'-^-' as defined by BS in [3^ as opposed to BL as was mentioned earlier . From these relations 
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and the explicit expressions provided for /i(i)'(3),(6) it is straightforward to evaluate the fields /i*-*^ and their 

inner and outer r-derivatives at r = tq, which then give us the total contribution of the monopole (£ = 0) to the GSF. 



4. The Even Static Modes (£ > 2(even), m = 0) 



These modes require a special discussion not only because the dimension of the homogeneous solutions space is 
smaller but also because the BC require extra care. With m = 0, we have that a*-^^ = and a*-*^ = . Furthermore, 
an inspection of tr,t6,t(j) components of ft-^^ (cf. Eq.(20) of [20|) reveals that these depend only on /i^^^ and h^^\ 
Since static modes must be symmetric under time reversal, we have that hu = for i = r,9,(j) thus we must have 
= and = for the static, even modes. This reduces the total number of fields in the even sector to 5 and 
eliminates the gauge equation ((T5|) (it gives the trivial = 0). Using the remaining two gauge equations ((T6| . (fT7| . 
we can obtain expressions for in terms of i?(i)'(3),(5)^ then substitute these into the field equations ([50)) - 

([32)) . This yields modified field equations for (i) = 1,3, 5: 



R 



(1) _ ;^(5) 



fR 



(3) 





f 


0- 






= Win, 








2/ f 






d^R^'^ 










+ 







/ 

47\f 



(50) 
(51) 

(52) 



Next, we calculate the boundary conditions for the static homogeneous solutions iij^g 5. Because we are looking at 
static modes, the ingoing/outgoing wave conditions are no longer appropriate for the BC. Our determining criterion 
is now regularity, so for the inner homogeneous solutions we select the following ansatz: 



R- 



CXD 

E 

— ^start 



(53) 



Substituting the ansatz ([53)) into the field equations for (i) = 1, 3, 5 gives us new recursion relations for the BC, which 
we display explicitly below as we will be making remarks about them here. We also list them in appendix lA 31 



8M^k{k - 2)bl 



^1-2^1-2 



'^k~2"k-2 

bl-3, 



2Mbl_^ 



+ Dl_A^^ + 2L{bl_.,-bl^^), 



f^3 l3 



-4Mbl_, 



SM^kbl 

- Fl_3bl_^ 



AM^bl_, 



fe-3i 



Di-2bl-2 



(54) 



(55) 



where 



El 



Ck 



-4M2(fc + l), Cl^4AI^k{k-l) 
2M{k~2), Dl^2M{L + k{l~2k)) 
L + l-k^, El = k-1, Gl^2Mk, 
AM'^{L + 1 + Ak- 3k^), Gl = 2Af (2L + 2 + 2fc - 3k'^). 



2ML-AM{l + k^) 
Dl^L-k{k + l), 



(56) 
(57) 
(58) 
(59) 



Little care is needed when evaluating the coefficients b\, b\ using the recursion relations ([54)) and ()55p . First, because 
the left-hand-side of Ea.([M)) gives zero for A; = 0, 2 we must start this recursion relation at /c = 3 with 6j = 6j = 0. 
Similarly, the recursion relation for b\ starts at fc = 2 with 69 = and b\ as the free parameter. Further inspection 
reveals that the remaining two free parameters are feg, 6?. This can be seen by realizing that Cj^ = for /c = 1 so we 
can not use the recursion relation ([55)) until k = 2 but we need 69 and bf to determine 6^>2 Eg. ([54)) and 5|>2 i'^ 
Eg. ([55)) . So our 3-dimensional basis of homogeneous solutions is generated from the set {fog, 6^, 6f }. When we evaluate 
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these three recursion relations to obtain the higher-k coefficients, we first get then at the (fc + 1)*^ order we 

recover 6|. For example, at fc = 2 we obtain 621^2 then at /c = 3 we recover 62 a-nd also obtain b^jb^. As usual, we 
truncate the infinite sum at some k = kmax such that the contribution of {kmax + l)*'^ term has absolute magnitude 
less than 10"". 

Next, we turn to determining the outer boundary conditions. This particular case is more involved than all the 
other BC thus far mentioned. First of all, the naive ansatz of = ^ka\/r^ only provides two free parameters 
thus falls one short of the needed three d.o.f. for the outer solutions. Inspired by the analytic, outer homogeneous 
solutions for I = odd, m = modes, which have r~^, Inr large-r behavior, we make the following ansatz 

^^f^aj^^ (60) 



When we substitute this ansatz into the ODEs (|50l). ([STj). (f52|) . wc find that a\ = for all k < £. Two of the three 
free parameters are a^^^^a^^^^ which combine to give 



a 



5 



^l = <^' + fh- (61) 



The next order terms in the recursion relations are as follows 

1 



[2L{2 + €)a\ - \La\ + 2(2 - e)a^f\ 



1^+1 = ^ {2L{1 - 2)a\ + 12La^ + + 2) - 2)af] , 



Note that all of these still only depend on the 2 free parameters a|, a|. It turns out the third free parameter is a|_|_2. 
As for the a^, they are all given in terms of {a^, a|, a^_^_2} with the condition a^^f_|_2 = 0. Unlike the previous cases, 
here we get two sets of recursion relations from each field equation, one for and another for a\.. These are: 

Clal = {k + l)al + al ~ 2kal ~ al 

-2M (Dl_,al_, + Dl_,al_, + + El_,dl-i - '^'4-1 



[F^_A-2 - 4-2) , (62) 

cl4 = {k + 1)4 + 4- 2M(^Di_,4_, + Dt,4^, + 4_,) 

+AM'Fi_24_„ (63) 



where 



Cl = L + l-k\ Dl^k{k-1), bl = 2{k + l), 
El ^ 1- 2k, = k + l. 



Clal = {k + l)al-al-4~2k4 

-2M {gI_,4_, + Gl_,4_, + HLi4-i - 24-i) 

+4Af2 (/t2aL2 + JL24-2) , (64) 
Cl4 = {k + l)4-4~2M{Gl_,4_,+Gl_,4_,) 

+4MHI_,4_,, (65) 

where 

Gl = 2^-2- L, Gl = 2k, Hi = -Ak, 
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Clal ^ 2L{al - al) - Dial + 2M (i^ti^Li + 2Lati + D^^al^,) , (66) 
Clal = 2L{al - 4) + 2M [eI_,4-i + 2L4_,) , (67) 

where 

Cl=L + k{l-k), bl = 2k^l, El^k{l-k) + 2. 

The careful reader will note that the recursion relations appear coupled to each other in Eqs. ([S^ - (|^7)) . That is, 
unlike all other recursion relations, the right-hand-sides of Eqs. - (p7)) contain k^^ order terms. If we move all 
order k terms to the left- hand- sides of Eqs. (|62|) - (|67l). we find that the LHSs form a coupled system of 6 equations 
with 6 unknowns. These equations are 'uncoupled' by using standard linear algebra methods. This naturally leads to 
the RHSs transforming into rather cumbersome expressions so we omit displaying them here. 

With the boundary conditions for the inner and outer homogeneous solutions computed, wc numerically solve 
the coupled set of three ODEs as before. The vector space of linearly independent homogeneous solutions is now 
6-dimensional and is constructed from inner, outer homogeneous solutions generated using BC obtained from the 
sets {6Q,6f,6f} for the inner and {a|, a^, a£_,_2} for the outer solutions, respectively. So at each (^ > 2,m = 0) even 
mode, we numerically integrate the ODEs for a total of 2 x 6 = 12 times. To determine the inhomogeneous solutions 
/j(i)X3),(5)j'j,^'j g^j^j their inner/outer r-derivatives, we construct a 6 x 6 complex matrix and invert it to solve for the 
complex amplitudes xi, . . . ,xq as before. We omit the details here as we have illustrated how to do this for both the 
generic odd and even modes in sections IIII A 11 IIIIB II respectively. Once these fields are known, wc can then use the 
gauge equations to construct i?(^^(ro) and R^'^^tq) and their inner/outer r-derivatives at r = tq. 



IV. COMPUTING THE GRAVITATIONAL SELF-FORCE 

With all the metric fields h!-"^^ and their r- derivatives computed, we now focus on the actual calculation of the 
gravitational self- force. We follow the prescription of [l| and [32j . 

Because we are modeling the small mass /i as a point particle, we are faced with the issue of the divergence of the 
GSF at the location of the particle. This requires a careful regularization of the GSF to remove the divergent, but 
non-physical, piece from it. We can write the regularized GSF as [l^ 

F'^ixo) = hm [FU{x) - F^U^)] , (68) 

X^Xq 

where Fj^jj is the "full" GSF constructed from the metric perturbation, and F^^^ is the "direct" (divergent) piece of 
it. Physically speaking, F^^^ can be thought of as representing the instantaneous part of the GSF that propagates 
along the past light-cone of the particle. 

In the mode-sum scheme, Ff^^ and F^^^ are decomposed into multipole modes Ff^/j and F^J. Thanks to this 
multipole expansion, the individual Z-modes of the divergent piece F^^j^. all have finite values at the x'^ — ?► Xq limit. I 
here represents the scalar spherical harmonic modes and it should not be confused with the tensorial modal index £ 
of the previous sections. 

Individual Z-modes of Ff^jj are obtained from the fields and their derivatives as given by Eg. ([7^ below. Then 

the GSF at the location of the particle (xq) is given by 

oo oo 

F"(xo) = ^ ([F£/,(xo)]± - - B-) ^ 5^[F,",^(xo)]± , (69) 

;=o /=o 

where L1/2 = 1 + 1/2. The ± correspond to taking the r-derivative at the r ^ Tq limit. A'^ and are regularization 
parameters. They are derived from the local structure of F^J near Xq^. + B" represents the asymptotic 

form of F^J for large I. For circular orbits in Schwarzschild spacetime, A^. = B" = for a = t, 9, (f). The non-zero 
r-components are given by [l3| 



^0 V ro J 



= ^ZJ[^- — ] ' (70) 
/i^roFg 



B 



E{w)-2K{w) , (71) 
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where Lq = (Afro)^/^/(l — 3Af/ro)^/^ is the orbital angular momentum, K{w) = /J^^^(l — wsin^ x) ^^'^dx and E{w) = 
/J^^^(l — w sin^ a;)^/^dx are the complete elliptic integrals of first and second kind, respectively and w = (ro/M— 2)^^. 
The regularized GSF can be computed by using either one of the ± values: the quantity F^^^"^ — Li/2^± is direction 
independent. This ± equality provides us with a way to check our GSF results. Since the i-component needs no 
rcgularization, we can write F^J^ ^ = j_ = F^J^^ _ . 
The I modes of the full force are given by [111 

[Kni^o)]^ - 7J E >^'"('^/2,'^o)x (72) 

m=-l 

[-7-a / — 3,m , ■j~al — 2.m. . -r-a l—l,m . -r-al.m . -r-a l-\-l.m . -r-a l+2,m . -r-a l+3,m 
•^(-3) +•^(-2) ' "'"■^(-1) "'"■^(0)' "'"■^(+1) "'"■^(+2) ' +-^(+3) 

■^S)™ ^^"^ constructed from /i^')^™, h!^]}^^ , Ji^f^"^ at x^ = Xq. The expressions for -T^^'™ are quite lengthy and are 
explicitly given in appendix C of [32| for circular and in appendix C of [l| for eccentric orlDits in Schwarzschild geometry. 
For this reason, we omit presenting them here. However, we would like to remark that •^(")™ contain coupling terms 

between tensor modes £ and scalar modes /. This is because the metric perturbation h^^ is decomposed in terms of 
tensor modes £, but the GSF is computed by summing over scalar modes / (the rcgularization procedure requires the 
mode decomposition to be done in spherical harmonics [l5| . [l8|). As a result, a given scalar spherical harmonic mode 
I will couple to 5 tensor spherical harmonic modes with £ — 2<Z<£ + 2for the r-component, and to 7 tensor modes 
£ — 3 < ^ < £ + 3 for the t-component of ^'^jy This is the reason why the index {j) in Eg. ([7^ goes from (—3) to (3). 

An extra simplification arises in Eg. ([7^ because the spherical harmonics F'™(7r/2, (/)o) = for I — m ~ odd. 
Furthermore, because F'™ — )- under m — > —to, we compute the sum only from m ~ \ to ra = (. then 

fold over the m-sum properly to include the m < contribution and finally add to these the to = term in the 
summation in Eg. (|72p . This is then regularized at each I mode via Eq. (|69p . 

To obtain the final value for the GSF, we compute the sum over all scalar I modes. Since the i-componcnt converges 
exponentially, Zmax ~ 10 suffices to obtain the value of F*{xo) to machine accuracy. However, the r-componcnt of 
the GSF falls off as L^^^ and this converges much more slowly. As we are using finite computer power to calculate 
an infinite sum over I, we must truncate the sum for the r-component at some / = Zmax (usually somewhere between 
15 and 30) and use fitting methods to estimate contribution from the I > Zmax modes. This contribution accounts 



for at most ~ 2% to the overall GSF [32| and is called "the large-Z tail" . The details of how to compute it arc given 
extensively in section HIE of [32| . Basically, one extrapolates the I > ^max terms in the sum using polynomial fits in 
powers of L^^^- As we use the same fitting method as [s^l, we refrain from elaborating any further. The details can 
be found there but let us discuss briefly how the tail error depends on the parameters used to do the fit. 

There are two free parameters that determine the large-^ tail. The first one is the number k of I modes G [/max + 
1 — A:, /max] that we select for the extrapolation. The second is A'^, which determines the degree of the polynomial fit 
m powers 

ofL-^^2- We use a numerical scheme that varies these two parameters {k, N) and finds the optimal values 
for both by comparing the error between the regularized I modes _F'^'^^n»''x+i-fc<i<imax obtained from the fitting formula 
and the actual numerical values computed by solving the Einstein equations. Our scheme uses the following ranges 
for the two parameters: 2 < A < 6 and 5 < fc < 12 depending on the total number of I modes that we compute 
(varies from 15 to 30). Because our frequency-domain code is able to compute up to 30 modes within an hour for 
tq < 20M, we are able to reduce the fractional error in the tail computation to ^ 10^®. As we will see below, the 
uncertainty in the largc-Z tail is not always the source of the most significant error in our computation. 



A. Summary of Methods and Computational Details 



Working in the frequency domain, we started by numerically solving the 10 coupled field equations for the 

radial fields i?^^(r) (the modes (£, to) = (0, 0), (odd, 0) have analytic solutions). To this end, for the first time, we 
calculated the boundary conditions for the radial fields in Lorenz gauge. We constructed linearly independent bases of 
homogeneous solutions and used these to obtain the inhomogencous solutions -^^^(ro) and their r r-derivatives 
via junction conditions. Following the prescription of [l^], we computed the •^(j)™^ i -^q)™^ of Eq. ([7^ . The I modes 
of the 'full' GSF are then given by this equation. We regularized the GSF at each I mode with the help of Eq. ((69|) 
then added all the individual Z-mode contributions together. Finally, for the r-component, we added the large-Z tail 
to the I sum to account for the F^J^'"'^'^'^ terms that we did not actually compute. It is this final result that equals 
the actual gravitational self-force. It is this quantity that we compare with BS in section fVl 
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Our numerical code is written in C and uses Gnu Scientific Library (GSL) repositories |52| for tfie numerical 
integration of the ODEs and matrix algebra used in obtaining the inhomogeneous solutions. After exhaustive numerical 
experimentation, we selected to work with the Runge-Kutta Prince-Dormand (rkSpd) numerical integration routine 
as this proved to be the fastest. For our matrix inversion, we opted for the lower- upper (LU) triangular matrix 
decomposition. We use a single desktop machine with two quad-cores to run our code, which proved to be more than 
sufficient for GSF computations for circular orbits. More than 95% of the computing time is taken up by the numerical 
integration of the coupled ODEs. This task is further multiplied because of the need to construct iV-dimensional bases 
of homogeneous solutions. For example, a GSF computation due to the first 15 scalar modes (i.e. tensor £ = . . . 18) 
numerically integrates various coupled ODEs a total of 2192 times. 

The speed of the numerical ODE integrator depends on a few freely specifiable parameters: the size of the integration 
domain ['''*„, r*„(], and the numerical accuracy thresholds (Arei,Aabs) used by the integrator. Given an ODE, the 
code picks the smaller of the two thresholds to integrate. We have empirically determined that a relative ODE solver 
accuracy of Aj-ei = 10"^" is sufficient for computing the GSF to within an overall fractional error of < 10"^ for runs 
with orbital radii 6M < ro < 50A/. However, for ro > 50M runs, we observed that Arci needs to be brought as close 
to machine accuracy as reliably possible i.e. 10~^^. This is because the transition region between the outer wave-zone 
(where the homogeneous fields ft,'^'' — > e~'"" and the region where the fields exhibit power-law growth (near 
To) is farther out for larger rp. Therefore, the numerical solutions can possibly grow by more than 20 orders of 
magnitude as the routine integrates from rout to tq. This fundamentally limits the accuracy that we can reach with 
a numerical integrator using double fioating point precision. After some numerical experimentation, we settled on 
a scheme that adaptively varies Arei, Aabs with increasing rg. The scheme works well for up to ro = lOOAf beyond 
which the accuracy thresholds thread very close to machine accuracy and the runtimes grow unreasonably long. 

The runtimes are rather insensitive to the location of r*^. The reason is that the potential Vim is very 'flat' near 
the event horizon (less than 1% variance as one goes from r*„ = —35 to —55), so the solutions hardly change. On the 
other hand, the runtimes do depend heavily on the location of r*^^. Therefore, its location must be chosen carefully. 
We elaborate more on this in the next subsection. 



B. The Error Budget 

The major sources of error that go into our computation are: (1) Error in the large-^ tail, (2) Error in the numerical 
matrix inversions used to construct the inhomogeneous solutions, (3) Numerical discretization error in the numerical 
integration of the ODEs, and (4) The fact that the boundary conditions are not computed at r* = ±oo. 

We determined that the error coming from the finitcness of the locations of the boundary points is much smaller 
than the other three sources of error. After some numerical experimentation, we came up with a satisfactory location 
for Tout (''out) keeping in mind the wave-zone condition rout >> (■To/tOm and the fact that our code slows down too 
much if Tout is unnecessarily too far out. This optimal choice was mentioned earlier in section Fill Al We tested the 
sensitivity of our solutions against changing Tout- We found that the relative variation in was < O(10~^^) when 
Tout was increased by up to one order of magnitude. 

We have already commented on the errors in the large-/ tail computation. Our usual standard has been a fractional 
error of 10~^ in the large-/ tail. As mentioned in section llVl we can reduce this error to nearly 1.0 x 10~^ by computing 
more numerical modes, but this naturally increases the runtimes. On the other hand, if we adhere to a fractional 
error of IQ-^ or 10"^ then we can reduce the overall runtimes considerably by computing less modes. We show this 
in Fig. m where we display plots of runtimes vs. rg for overall fractional errors of 10~^, 10~^ and 10~^. In short, we 
have a good understanding and good control over the uncertainty in the large-/ tail. 

The numerical discretization error coming from the numerical integration of the ODEs contributes much less to 
the overall error than the other error sources mentioned here. The GSL ODE integrator routines are very robust and 
have a very good handle on discretization errors. Our own numerical tests showed that these errors have magnitudes 
< 0(10"^^) with respect to the inhomogeneous fields. 

Finally, as mentioned in section flllB 11 the biggest source of error comes from the numerical inversion of the matrix 
constructed from the homogeneous solutions. This becomes the dominant source of error for rg > 50Af . An inspection 
of the matrix inversion output for each (£, m) mode reveals that the inversion errors grow with increasing £ — m and 
that they are also larger in the even parity sector. We monitored the condition numbers of the matrices and found 
out that for even parity modes with £ — m > 15, they routinely exceeded 10^^ for ro > 50A/ and got as large as 
10^^ for To > lOOA'/. Further inspection of these large £, large rg even modes revealed that the determinant threads 
very close to zero. This is an indication that our linearly independent bases of homogeneous solutions start becoming 
degenerate in this region. The reason why this happens for large £ — m is due to particular way we have formulated 
the location of the outer boundary by setting rout = 50fro/a;m = 50rQ^'^{£/m). From this, one sees that rout reaches 
its maximum value when £ — m reaches its maximum value. So, this 'degeneracy problem' is actually caused by large 
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values for rout • What happens is that because the leading order power-law for each homogeneous field dominates near 
To, the solutions that have the same power-law behavior start looking numerically identical as the integrator works 
its way in toward tq. As we look at the values of the fields for larger rg runs, the matrices constructed from the even 
parity homogeneous fields become linearly dependent (singular valued). This means the matrix inversion is not very 
reliable. We find that this degeneracy of even parity solutions becomes significant for the runs where tq > 50Af . So, 
any numerical ODE integration that routinely goes beyond this point (ro ~ 50M) starts running into this degeneracy 
problem. 

We model the error coming from the singular-valucdness of the matrices as a continuity violation in the inhomo- 
geneous fields i?(*)(r) at r = ro . This continuity violation, A^*), is most prominent for the fields i?(2),(4),(5) -^yjiCTc it 
is about 0(10^) larger than the violations in the other fields. In the worst case, e.g. ro = 150M and ^ = 17, m = 1; 
A(5) « 10^^ However, even at rg = 150M, the violation quickly subsides to < 10 ^ once m > 2 whatever £ may be, 
but because the GSF is constructed by summing over all {£, m) modes, this error is additive. For a computation of 
the GSF requiring ^max = 18, the relative strength of the error is amplified by a factor of ~ 10^ — 10'^ going from from 
a single mode to the final GSF, which is constructed from the sum of 0(10^) modes. This is indeed what we observe 
numerically. We have not yet looked into fixing this inversion problem but we are aware that using singular- valued 
decompositions for the matrices do not offer an improvement |53| . Be that as it may, we do not think this to be a 
problem for when we compute the GSF for eccentric orbits because we will be mostly interested in the strong field 
regime of ro < 20M. However, for equatorial eccentric orbits in frequency domain, we expect to encounter a similar 
type of degeneracy in our solutions due to the fact that the frequency spectrum is determined by two fundamental 
frequencies: LUmn = rnil^ ~f nflr- There will be points in the parameter space where the two terms in LUmn will 
conspire to cancel each other to values less than 10~^. When this happens, the conditions numbers for matrices of 
homogeneous solutions grow to values that render the matrix inversion unreliable. We are currently working on a 
solution to this problem. 



V. RESULTS 



We present the output of our frequency-domain code for the gravitational self- force in Tables 2 and 3. For com- 



parison, we include the results of BS [32[ and the relative difference between our respective values for the t-, and 
r-components of the GSF. We find very good agreement with the results of BS (within their error bars) for ro up 
to ~ lOOA/. However, beyond that, our values stray from theirs. Given that Berndtson agrees with BS within 
their quoted errors bars for up to 150M, we must conclude that the degeneracy problem renders our results unreliable 
beyond ro ^ lOOAf . However, as our results show, in the strong field regime our f-domain results are much more 
accurate than their time-domain counterparts. 

As another way of confirming our results and determining the magnitude of the error in our GSF computation, we 
computed the energy flux of the gravitational waves leaving the system and compared the value of the total radiated 
power with the total rate of energy loss given by the dissipativc component of the GSF. In the case of circular orbits, 
only the t-component of the GSF is dissipative so the rate of energy loss can be related to i^* as follows: 

— = Ft. (73) 

In terms of Schwarzschild time t, this becomes dEo/dt = — (/xuo)~^-Fi, where Uq is the t-component of the 4-velocity 
of the particle evaluated at r = rg. In the adiabatic approximation, where fJ-/M 1, dEg/dt can be taken to be the 
average rate of energy loss per orbit. Energy conservation dictates that this loss of energy must be balanced by the 
total energy flux carried by gravitational waves radiated out to infinity and absorbed into the black hole. Therefore, 
we have the following balance equation: 

^totai = E^ + Eeh = = Ft/ul, (74) 

dt 

where the overdot now denotes d/ dt and Eoo, ^-eh denote the gravitational wave flux radiated to infinity and through 
the event horizon (EH), respectively. These fiuxes are constructed from the metric fields ft,'^*^^™. We omit the details 
of this construction here, but for the interested reader they can be found in [H, [13, HI] . Let us simply display the 
final expressions for the fluxes: 



°° 647rA£(£-hl) 

£=2 rn=-l ^ ' 



(75) 
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2567rM2(l + 16M2m2rj2) 



-(1) l + 4iMmrJo 



where ft.|^^EH implies that the metric field is to be evaluated (in frequency domain) at r = cxj, 2Af , respectively. Using 

Eqs. (|75p . (j76p and our results for i^*, we compute the total radiated power i^totai and compare the resulting values. 
The relative difference between the two results is shown in the last column of Table 3, which shows that the agreement 
is excellent for small tq. It naturally gets worse for increasing values of tq. We also found that the disagreement 
between the two values for iJtotai matched our overall fractional error in well. 

We also present the runtimes for our code for three different relative accuracies. These are quantified by the overall 
fractional error in our numerical computation of the GSF. We have selected to present results for overall fractional 
errors of 10~*, 10~^, 10~^. We display the runtimes for these in Fig. [T] As can be seen from the upper left panel 
of the figure, at a relative accuracy of 10~^, our code takes less than two minutes to compute the GSF for radii less 
than ^ 15M. This grows nearly to a day as tq approaches lOOAf . Although toward lOOAf the runtimes appear to 
level off, this is due to our logarithmic scale for the vertical axis. The runtimes increase by ^ 100 minutes in going 
from 70A/ to 80Af, and 80Af to 90Af . In the same figure, upper right panel, one sees that demanding an accuracy 
of 10^^ increases the runtimes by a factor of two to three for ro < lOAf. However, beyond tq = 50Af, this accuracy 
becomes unattainable. Finally, we find it quite difficult to keep the overall fractional error less than 10""^. But as 
the lower left panel of the figure shows, an accuracy standard of 10~^ is achievable for tq < 30A^ and the overall 
runtimes are not prolonged by much for these strong field GSF computations. Interestingly enough, in the regime 
''o ^ 20Af, the To < 8Af runs seem to take more time than ro > 9Af runs. This was artificially caused by our need 
to compute more modes in order to lower the error in the large-/ tail for the tq < 8A/ runs. It turns out that for 
the smallest radn, the large- Z tail can not be computed to the desired accuracy of 10~^ or 10"'' using just 15 or 17 
scalar modes, which is what we had done for the ro > 9A^ runs. We think the reason for this is that the magnitudes 
of the individual I modes of the GSF are large enough for ro < 8Af that more modes are needed in order for the 
tail to be fit correctly. Finally, in the lower right panel, we present the computation times for a given ro < 20Af 
run for all three accuracies. As expected, the runtimes increase with demand for higher accuracy. However, by how 
much they increase is not the same at each radius. There is also the anomalous data point for the lOAf run where 
the 10~^ accuracy computation takes slightly less time than the 10~^ one. This comes from our not having explored 
thoroughly enough the free parameters that determine the overall error and runtime such as /max, number of points 
used in the tail and the numerical ODE integrator accuracy thresholds. Most importantly, the figure shows that all 
ro < 20Af runs take less than 15 minutes up to an accuracy of 10~^. 

It should also be added that even on our modest desktop, we can simultaneously perform a dozen strong field runs 
without significantly affecting individual runtimes. For example, in a 15 minute period, we can compute the GSF for 
all integer orbital radii from 6A/ to lOAf to an accuracy of 10~^. We find the speed of our code to be fast enough 
to encourage continuing this frequency-domain approach to tackle the eccentric Schwarzschild problem for the GSF. 
Work is currently underway and the preliminary results are encouraging. We intend to apply these methods to the 
full Kerr problem later on. 
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Appendix A: The Recursion Relations for The Boundary Conditions 

1. Generic Odd and Even Modes 

Recall that by generic, we mean the non-static (m 7^ 0), € > 1 modes. Here, lo denotes = mfio a-nd L = £{£ + 1). 
We begin by redisplaying the recursion relations for the outer boundary conditions (BC) for odd parity homogeneous 



20 



ro/M 




(M//i)2F5s 


Rel. 


diff. 


6.0 


2.4466495(4) x 10^^ 


2 


44661 


X 


10" 


2 


4 





X 


10" 


6 


7.0 


2.149907776(8) x 10"^ 


2 


14989 


X 


10" 


2 


8 


3 


X 


10" 


6 


8.0 


1.8357830(4) x 10"^ 


1 


83577 


X 


10" 


2 


7 


1 


X 


10" 


6 


9.0 


1.5637099(1) X 10"^ 


1 


56369 


X 


10" 


2 


1 


3 


X 


10" 


5 


10.0 


1.3389470(2) x 10"^ 


1 


33895 


X 


10" 


2 


2 


2 


X 


10" 


6 


11.0 


1.155174593(6) x 10"^ 


1 


15518 


X 


10" 


2 


4 


7 


X 


10" 


6 


12.0 


1.00462381(8) x 10"^ 


1 


00463 


X 


10" 


2 


6 


2 


X 


10" 


6 


13.0 


8.8048853(3) x 10"^ 


8 


80489 


X 


10" 


3 


5 


3 


X 


10" 


7 


14.0 


7.7730602(4) x lO"'' 


7 


77307 


X 


10" 


3 


1 


3 


X 


10" 


6 


15.0 


6.9081719(3) X 10"^ 


6 


90815 


X 


10" 


3 


9 


7 


X 


10" 


5 


20.0 


4.1570550(2) x 10"^ 


4 


15706 


X 


10" 


3 


1 


2 


X 


10" 


6 


30.0 


1.9698169(3) x 10"^ 


1 


96982 


X 


10" 


3 


1 


6 


X 


10" 


6 


40.0 


1.142883(1) x 10"^ 


1 


14288 


X 


10" 


3 


2 


6 


X 


10" 


6 


50.0 


7.449480(1) x 10~* 


7 


44949 


X 


10" 


4 


1 


3 


X 


10" 


6 


60.0 


5.236083(3) x 10"* 


5 


23613 


X 


10" 


4 


9 





X 


10" 


6 


70.0 


3.8801(1) X 10"* 


3 


88010 


X 


10" 


4 


2 


6 


X 


10" 


6 


80.0 


2.9896(1) X 10-* 


2 


98979 


X 


10" 


4 


6 


4 


X 


10" 


5 


90.0 


2.3739(1) X 10"* 


2 


37406 


X 


10" 


4 


6 


7 


X 


10" 


5 


100.0 


1.9304(1) X 10~* 


1 


93063 


X 


10" 


4 


1 


2 


X 


10" 


4 


120.0 


1.3483(1) X 10"* 


1 


34868 


X 


10" 


4 


2 


8 


X 


10" 


4 


150.0 


8.673(1) X 10"^ 


8 


68274 


X 


10" 


5 


1 


1 


X 


10" 


3 



TABLE II: Output for the r-component of the gravitational self-force for various orbital radii ro compared with results of BS 
[3^ . Column 2 contains our results; the number in parentheses indicates the size of the uncertainty in the last significant 
digit, e.g. 2.4466495(4) = 2.4466495 ± 4 x 10"^. In column 3, we display the resuhs of BS for comparison. Column 4 gives 
the relative difference between our values and BS'. Our results are within their quoted error bars for nearly up to ro = lOOM. 
Beyond that the disagreement seems to grow up 0(10"''). Given that Berndtson's results [35| agree with BS better for large 
ro, we conclude that our current results are not reliable beyond ro ~ lOOM. Nevertheless, as can be seen from the number of 
significant digits that we have included for for ro < 50M, the frequency-domain results are much more accurate than time 
domain in the strong field regime. 



fields Rt and R^^. 



2i(jjk at 



2iujk al" 



Ck-i al_, + Dk-2 4-2 + Ek-3 4-3 + 24-1 - lOMfli'i^ + UM'al^L,, 
h-i + Jk-2 alH^ + Kk-3 all, + 2Xal_, - UIXal_^, 



where 



Ck = Uliujk + fc(fc + 1) - L - 4, Ik= Uliojk + k{k + 1) - L + 2, 

Dk = -6Mk - AMk^ + 247\/ + 2ML, Jk = ~6Mk - Ulk^ - 6M + 2ML, 

Ek = 4A/2(/s2 + 2/c-8), A'fc = 4M2(/fc2-|_2fc + l). 

Next, we present the recursion relations for the inner BC 

AM^k{k~AMi^)hl ^ Ck-ibl_^+Dk-2bl_2 + Ek-3 + 2Mhf_^ ~ 2b]^_2, 
4M^k{k ~ AMiuj)bl'' = Hk-i + Jk-2 fofc°_2 + bll, - U4\bl_^ - 2A5; 



fc-2 



where 



Ck = ^Mik + 12Miujk -2k^ + L- 4), Hk = 2M{k + 12Miujk - 2k^ + L - 1), 
Dk = 4+ UMiujk + L- k{k - 1), Jk = -2 + UMiuk + L- k{k - 1), 
Ek — 2iujk. 



(Al) 
(A2) 



(A3) 
(A4) 
(A5) 
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ro/M 


(M//i)2F* 




t 

BS 




Rel. 


diff. 


{M/^lfEtotal 


{M/iifF,/ui 


Rel. diff. 


6.0 


-1.9947610064(3) x 10^^ 


-1.99476 


X 


m 


-3 


5.0 X 


10' 


7 


9.4033935631 x lO"* 


9.4033935626 x 10"* 


5.7 X 10" 


10 


7.0 


-7.411127850(9) x 10""* 


-7.41101 


X 


10- 


-4 


1.2 X 


10" 


5 


4.001632906 x 10"" 


4.001632909 x 10"* 


6.6 X 10" 


11 


8.0 


-3.307397510(3) x lO"* 


-3.30740 


X 


10" 


-4 


1.2 X 


10" 


5 


1.9610454858 x 10~* 


1.9610454864 x 10"* 


3.0 X 10" 


10 


9.0 


-1.668101230(4) x 10"* 


-1.66810 


X 


10" 


-4 


1.2 X 


10" 


5 


1.0593325177 x 10"* 


1.0593325178 x 10"* 


8.8 


X 10" 


11 


10.0 


-9.19075772(7) x 10"^ 


-9.19067 


X 


10- 


- 5 


1.0 X 


10" 


5 


6.151631678 x 10"^ 


6.151631677 x lO"'^ 


2.2 X 10" 


10 


11.0 


-5.41623002(6) x 10~^ 


-5.41605 


X 


10" 


-5 


1.2 X 


10" 


5 


3.779162580 x 10"^ 


3.771962578 x 10"^ 


4.8 X 10" 


10 


12.0 


-3.3659568(1) x 10"^ 


-3.36587 


X 


10" 


-5 


1.2 X 


10" 


5 


2.42917009 x 10"^ 


2.42917010 X 10"^ 


3.2 


X 10" 


-9 


13.0 


-2.1839249(2) x 10"^ 


-2.18388 


X 


10" 


-5 


1.2 X 


10" 


5 


1.620747493 x 10"^ 


1.620747489 x 10"^ 


2.8 


X 10" 


-9 


14.0 


-1.4685410(2) X 10"^ 


-1.46851 


X 


10" 


-5 


1.2 X 


10" 


5 


1.115762106 X lO"-"^ 


1.115762104 X 10"^ 


2.2 


X 10" 


-9 


15.0 


-1.0177145(1) X 10"^ 


-1.01772 


X 


10" 


-4 


1.2 X 


10" 


5 


7.88902019 x 10"*^ 


7.88902015 x lO"*^ 


5.3 


X 10" 


-9 


20.0 


-2.2554391(2) x lO"*^ 


-2.25549 


X 


10" 


-6 


2.2 X 


10" 


5 


1.87147091 X 10"^ 


1.87147088 X 10"*^ 


1.6 


X 10" 


-8 


30.0 


-2.8081894(8) x 10"^ 


-2.80813 


X 


10" 


-7 


2.1 X 


10" 


5 


2.486484 x 10"^ 


2.486486 x 10"^ 


1.1 


X 10" 


-6 


40.0 


-6.51228(2) X 10"** 


-6.51219 


X 


10" 


-8 


1.5 X 


10" 


5 


5.95014 X 10"* 


5.95015 X 10"* 


2.9 


X 10" 


-6 


50.0 


-2.108456(4) X 10~** 


-2.10849 


X 


10" 


-8 


3.0 X 


10" 


5 


1.962458 X 10"** 


1.962453 X 10"* 


2.3 


X 10" 


-6 


60.0 


-8.41300(3) X 10"^ 


-8.41306 


X 


10" 


-9 


7.0 X 


10" 


6 


7.92644 X 10"'-* 


7.92641 X 10"^ 


3.9 


X 10" 


-6 


70.0 


-3.8743(1) X 10"^ 


-3.87411 


X 


10" 


-9 


4.1 X 


10" 


5 


3.6819 X 10"^ 


3.6818 X 10"^ 


3.5 


X 10" 


-5 


80.0 


-1.9804(1) X 10"^ 


-1.98069 


X 


10" 


-9 


4.6 X 


10" 


5 


1.8945 X lO"'' 


1.8946 X 10"^ 


4.2 


X 10" 


-5 


90.0 


-1.0966(3) X 10"^ 


-1.09654 


X 


10" 


-9 


9.1 X 


10" 


5 


1.0541 X 10"^ 


1.0544 X 10"^ 


2.6 


X 10" 


-4 


100.0 


-6.464(2) X 10"^" 


-6.46305 


X 


10" 


-9 


2.1 X 


10" 


4 


6.238 X 10"^° 


6.240 X 10"*° 


3.2 


X 10" 


-4 


120.0 


-2.596(9) X 10""' 


-2.59096 X 10" 


10 


1.9 X 


10" 


3 


2.516 X 10"^" 


2.525 X 10"*° 


3.6 


X 10" 


-3 


150.0 


-8.44(6) X 10"" 


-8.47172 


X 


10" 


11 


3.4 X 


10" 


3 


8.27 X 10"" 


8.22 X 10"*^ 


6.1 


X 10" 


-3 



TABLE III: Output for the t-component of the gravitational self- force for various radii ro. Our results are in column 2. In 
column 3, we show the results of BS for and display the relative difference in column 4. Once again, our results fall within 
BS' error bars for ro < lOOM and again the error increases up to 0(10""^) for ro = 150M. And as was the case with F"^ , 
our frequency-domain results for F* also have much smaller uncertainties in the ro < 50M regime compared with those of BS. 
We also checked our results for using energy balance arguments. Since only the t-component of the GSF is dissipative for 
circular orbits, it can be related to the energy flux leaving the system as we have outlined in section |V] The total energy flux 
is computed using the two different methods and these results are displayed in columns 5 and 6 down to the significant digit 
at which they start disagreeing. Column 7 contains the relative difference between the two values. Once again, the agreement 
is extremely good for small ro and grows to 0(10"'') as ro increases to 150M. 



Now, we turn our attention to the BC for the even parity fields 5 6 7- start with the recursion relations for 
the outer BC for Rf,Rf and Rf: 

2iwkal = + (2 - 4Afi^)aLi + 2aLi + 2al_^ + Dl_^al_^ + Dl_^al_^ - 12Afat2 



where 



(A6) 

2iwkal = Cl^A^i + 2(4_i - - a^i) + DI-^^^^ 

+ 4A/(-aL2 + 4-2 + 3«L2) + El_,al_, - 16M^at_„ (A7) 
2tujkal = Cl_,al^, + 2(4_i - 4-1 - 4_,) + DI_^4_, 

+ AM i-4_^ + 4.2 + ^4-2) + - ^QM^4-3, (A8) 

d = k{k + 1)+ AMiujk - 2 - L, (A9) 

Dl = 2Af(5 + L-2fc2 -3fc), L*^. = 2A/(2fc - 8 + 4A/ia;), (AlO) 



'k — ^ivi \^ ^ 1^ - ^n, -OA,;, ijf^ 

El = AM^{k^ +2k-i), El = 8M^{5~2k), = 16M^{k - 2). (All) 



For the field Rf, we have: 



2tu:k4 = Cti4-i + 2^(4-1 - 4-1 - 4-i) + K-1 + Dl_A-2 - lOAfaLa (A12) 
+ 2ML(-24_2 + 4at2 + ^4-2) + EI_^4-z + AM^{-2L4_^ - SLatg + 34 
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Strong Field Runtimes vs. Accuracy 
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FIG. 1: The runtimes for the 10"'', lO"*^, 10"'^ overall fractional error runs. Panels (a), (b) and (c) display plots of runtime (in 
minutes) versus orbital radius ro at which we compute the GSF. Aqsf denotes the overall fractional error in our numerical 
computation of the GSF. This error is what we refer to as our (relative) 'accuracy'. As can be seen in panel (a), at an accuracy 
of 10~*, our code takes less than two minutes to compute the GSF for ro < 15M. This grows nearly to a day as ro approaches 
lOOM. Panel (b) shows that an accuracy of 10~® increases the runtimes by a factor of two to three for ro < lOM, but the 
runtimes are still < 10 minutes for ro < 20M. However, beyond ro = 50M, this accuracy becomes unattainable. As panel 
(c) shows, an accuracy of 10~^ is achievable for ro < 30M and the overall runtimes do not change much for these strong field 
GSF computations. Interestingly enough, for ro < 20M, the ro < 8M runs seem to take more time than ro > 9M runs. This 
is a result of our having to compute more modes to obtain the GSF for the ro < 8M runs because the large-Z tail could not 
be computed to the desired accuracy of 10~® or 10~^ using just 17 scalar modes, which is what we had done for the ro > 9M 
runs. We think the reason for this is that the magnitudes of the individual I modes of the GSF are large enough for ro < 8M 
that more modes are needed in order for the tail to be fit correctly. Finally in panel (d), we present the runtimes for a few 
ro < 20Af run for all three accuracies. As expected, the runtimes increase with demand for higher accuracy (except for the 
lOM run). Most importantly, the figure shows that all ro < 20Af runs take less than 15 minutes up to an accuracy of 10~^. 



where 



And for R^: 



Cl = k{k + l)+4Miujk-4:- L, Dl = 2M{12-2k'^ -3k + L), (A13) 
El = 4:M^{k^ +2k-8). (A14) 



2iu;kal = Cl^,al_,+2\4_^+Dl„^al_^-4MXal_^+El_3al_3, (A15) 

where 

Cl = k{k + l)+4:Miu}k- L + 2, Dl^2M{L-3-2k^ -3k), (A16) 

El = 4Af2(fc2 + 2fc + l). (A17) 
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Next, we present the recursion relations for the inner boundary conditions for the same fields in the same order. We 
start with: 

8Af3fc(4M*w - k)bl = + CLi&Li - SAf^^ti + + Dl-,bl^2 (A18) 

-MIbl_^ + El^A-s + 2(1 + 2Anu;)bl_, + 2bl„, + 2bl_, - ^,^461-4, 



where 



Cl = AM^il + 'Sk^ - L-16Miiuk-k), Cl = 8M^i2Miuj - k), 
Dl = 2A/(3fc2 - 2fc - 2L - 1 - iUIiujk), = AM{AMiuj -k-1), 
El = k{k - 1) - UMitok -2-L, = 2iuk. 



(A19) 
(A20) 
(A21) 



AM^iAAHu - k)bl = GLi5ti+4Af(feti-^Li+d) + ^L2&L2 



(A22) 



where 



3 ^,6 

fc-2 



4M^k{4A'ULo - k)bl = Gl_,bl_^ + AMibU - bl_, + bl_,) + Hl_^b\ 

Gl = 2M(2fc2 - fc - L + 1 - UMiuk), 
Hi = fc(fc - 1) - i - 2 - l2Miujk. 



(A23) 

(A24) 
(A25) 



where 



And, finally 



where 



4Af2fc(4Mzc. - k)bl = /tifeti + 2ML{2bl_, + bl_{) - 2Mbl^, + JI_A_^ 
+ 2L{bl_, - bl_, - bl_^) + 2bl_, - ^La&ts, 

il = 2M{2k^ -k- L + 4:-12Miujk), 
Jl = - f ) - L - 4 - UMiuk. 



4:M^k{AMiuj - k)bl = Gl_j^bl_j^+AMXbl_j^+Kl_2bl_2 + 2Xb 

~ ^fc-3^fc-3i 



KI^ k{k-l) - L + 2- 12Miu}k. 



5 

k-2 



(A26) 

(A27) 
(A28) 



(A29) 
(A30) 



2. The Even Dipole (£ = l,m = 1) Mode 



The recursion relations for aj.' ^J.''^'^ do not change. However, we end up with new recursion relations for the 
inner and outer boundary conditions for Rf 

4M^k{4Miu; - k)bl = iLibl-i + 2ML{2bl_, + bl_,) - 2Mbl_, + Jl_2bl_2 



(A31) 



2iujkal = Cf_ia|_i+2i:(4_i-a3_^_a6_j^2aI._i+i?,^_2-fc-2 
+ 2A'/L(-24_2 + 4aL2 + Sa^s) + ^L3«L3 
+ 4M2(-2Lat3-3Lat3 + 3aL3), 



fe-2 



(A32) 



where w = mfJo and L = t{l + 1) as before. The coefficients C^, _D^., E^,F^,I^, are the same as before, displayed 
in Eqs. (UT31) . (UTl) . (IMTj) . (IM7)) and (1X281) . 
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3. The Static (m = 0) Even Modes 

Recall that because we are now dealing with static modes, we no longer have in or outgoing waves at the boundaries. 
For the inner BC, we use the following ansatz: 

oo 

Ri^ E ^fc(^-2Mf, (A33) 

A; — fcstart 

where we now have three fields labelled by i = 1, 3, 5. New recursion relations for the inner BC are: 

M'I^k{k-2)bl = FLibl-i+Gl_,bl_, + Gl_,bl_,-2Mbl_^ 

+EL3bLs + El_,bl_, - bl_„ (A34) 
4Mk{k-l)bl = Cl_,bl_,-AMLbl_,+Dl_,bl_^ + 2L{bl_^~bl_^), 
Cl_,bl_, = Cl_,bl_, - SM^kbi + AMHI_, + Dl_,bl_, + Dl_,bl_^ 

+UIbl_^ + El-3bl-3 + EL3bi-3 + bL3, (A35) 

where 

Cl ^ -4M^{k + l), Cl^AM^k{k^l), Cl^2ML-AM{l + k^), (A36) 
Dl ^ 2]\I{k-2), Dl^2M{L + k{l-2k)), Dl^L-k{k + l), (A37) 
El = L + l-k^, El^k-1, Gl^ 2Mk, (A38) 
= 4:M^{L + l + 4:k-3k^), 01 ^ 2M{2L + 2 + 2k ~ 3k^). (A39) 

For the outer boundary conditions, we make the following ansatz: 

Recall that the recursion relations for a],,a)l,a^,a]^,af.,a^ are determined by the three free parameters af,a^,a^_^2- 
We now present these relations in their coupled form: 

Clal = {k + l)al + al - 2kal - al 

-2M (^LiaLi + i^Liati + ati + i^LioLi - 2aLi) 

+4A/2 (f3_2«L2 - , (A41) 

Cldl = {k + 1)4 + 4- 2M(^Dl_,al_,+Dl_,al_,+al 

+AM'Fi_24_,, (A42) 



where 



Cl = L+l-k\ bl = k{k-l), Dl = 2{k + 1), 
El = l-2fc, Fi^k + 1, 



OkO-k — 



(fc + 1)4 - 


-4- 


a\ — 2fc4 


-2M (gL 






+4A/2 (/3 


3 


+ ^-2«i 


(fc + 1)4 - 




2 A/ (gLi 


+4A/2/3_2 







where 



Clal = (fc + l)4-4-2A/ Gtiati+Gtial 



G^ = 2fc2 _ 2 - L, G^ = 2fc, Hi = -4fc, 
/3 = fc2 - 1, J3 = -2fc. 



(A43) 
(A44) 
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Clal ^ 2L{al - 4) - ^|a| + 2M [eLA^, + 2Lal^, + DtA-,) , (A45) 
Cl4 = 2L{al - al) + 2M (eI_A-i + 2LaLi) , (A46) 

where 

Cl = L + k{l-k), bl^2k-l, El^k{l-k) + 2. 



[1] L. Barack, N. Sago, Phys. Rev. D81, 084021 (2010), [arXiv: 1002.2386 [gr-qc]]. 

[2] J. L.Barton, D. J. Lazar, D. J. Kennefick, G. Khanna, L. M. Burko, Phys. Rev. D 78, 064042, (2008), [arXiv:0804.1075 
[astro-ph]]. 

[3] [http: / /www.hgo. caltech.edu /advLIGO /scripts /summary.shtml] . 
[4] [http: / /wwwcascina. virgo.infn.it / advirgo /] , 

[http://wwwcascina.virgo.infn.it/advirgo/docs/whitepaper.pdf]. 
[5] [http://hsa.gsfc.nasa.gOv/Docunientation/LISA-LPF-RP-0001_vl.l.pdf]. 
[6] [http://lisa.nasa.gov/]. 

[7] F. D. Ryan, Phys. Rev. D 56, 1845 (1997). 

[8] S. Finn, K. Thorne, Phys. Rev. D 57, 7089 (1998). 

[9] P. Amaro-Seoane, J. R. Gair, M. Freitag, M.C. Miller, L Mandel, C. J. Cutler, S. Babak, Class. Quant. Grav. 24, R113 
(2007) , [arXiv:astro-ph /0703495] . 
[10] L. Barack, C. Cutler, Phys. Rev. D 69, 082005 (2004), [arXiv:gr-qc/0310125]. 
[11] B. S. DeWitt and R. W. Brehme, Annals Phys. 9, 220 (1960). 

[12] Y. Mino, M. Sasaki and T. Tanaka, Phys. Rev. D 55, 3457 (1997), [arXiv:gr-qc/9606018]. 
[13] T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381 (1997), [arXiv:gr-qc/9610053]. 
[14] S. Detweiler and B. F. Whiting, Phys. Rev. D 67, 024025 (2003), [arXiv:gr-qc/0202086]. 
[15] L. Barack and A. Ori, Phys. Rev. D 61, 061502 (2000), [arXiv:gr-qc/9912010]. 
[16] L. Barack, Phys. Rev. D 64, 084021 (2001), [arXiv:gr-qc/0105040]. 

[17] L. Barack, Y. Mino, H. Nakano, A. Ori and M. Sasaki, Phys. Rev. Lett. 88, 091101 (2002), [arXiv:gr-qc/0111001]. 

[18] L. Barack and A. Ori, Phys. Rev. D 67, 024029 (2003), [arXiv:gr-qc/0209072]. 

[19] L. Barack and A. Ori, Phys. Rev. Lett. 90, 111101 (2003), [arXiv:gr-qc/0212103]. 

[20] L. Barack and C. O. Lousto, Phys. Rev. D 72, 104026 (2005), [arXiv:gr-qc/0510019]. 

[21] L. M. Burko, Class. Quant. Grav. 17, 227 (2000), [arXiv:gr-qc/9911042]. 

[22] L. M. Burko and Y. T. Liu, Phys. Rev. D 64, 024006 (2001), [arXiv:gr-qc/0103008]. 

[23] L. Barack and L. M. Burko, Phys. Rev. D 62, 084040 (2000), [arXiv:gr-qc/0007033]. 

[24] L. M. Burko, Phys. Rev. Lett. 84, 4529 (2000), [arXiv:gr-qc/0003074]. 

[25] S. Detweiler, E. Messaritaki and B. F. Whiting, Phys. Rev. D 67, 104016 (2003), [arXiv:gr-qc/0205079]. 
[26] L. M. Diaz- Rivera, E. Messaritaki, B. F. Whiting and S. Detweiler, Phys. Rev. D 70, 124018 (2004), [arXiv:gr-qc/0410011]. 
[27] W. Hikida, S. Jhingan, H. Nakano, N. Sago, M. Sasaki and T. Tanaka, Prog. Theor. Phys. 113, 283 (2005), [arXiv:gr- 
qc/0410115]. 

[28] R. Haas, Phys. Rev. D 75, 124011 (2007), [arXiv:gr-qc/0704.0797]. 

[29] I. Vega and S. Detweiler, Phys. Rev. D 77, 084008 (2008), [arXiv:0712.4405]. 

[30] L. Barack and C. O. Lousto, Phys. Rev. D 66, 061502 (2002), [arXiv:gr-qc/0205043]. 

[31] T. S. Keidl, J. L. Friedman and A. G. Wiseman, Phys. Rev. D 75, 124009 (2007), [arXiv:gr-qc/0611072]. 
[32] L. Barack and N. Sago, Phys. Rev. D 75, 064021 (2007), [arXiv:gr-qc/0701069]. 
[33] S. Detweiler, E. Poisson, Phys. Rev. D 69, 084019 (2004), [arXiv:gr-qc/0312010]. 
[34] S. Detweiler, Phys. Rev. D 77, 124026 (2008), [arXiv:0804.3529 [gr-qc]]. 
[35] M. V. Berndtson, Ph.D. Dissertation, [arXiv:0904.0033 [gr-qc]]. 

[36] N. Sago, L. Barack and S. Detweiler, Phys. Rev. D 78, 124024 (2008), [arXiv:0810.2530 [gr-qc]]. 
[37] E. Poisson, Phys. Rev. D 52, 5719 (1995); Phys. Rev. D 55, 7980 (1997). 
[38] N. Warburton, L. Barack (2011), [arXiv:1103.0287 [gr-qc]]. 

[39] L. Barack, A. Ori, N. Sago, Phys. Rev. D 78, 084021 (2008), [arXiv:0808.2315 [gr-qc]]. 
[40] E. Poisson, Living Reviews in Relativity 7, 6 (2004). 

[41] L. Barack, Class. Quan. Grav. 26, 213001 (2009), [arXiv:0908.1664 [gr-qc]]. 
[42] T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957). 
[43] V. Moncrief, Ann. of Phys. 88, 323, (1974). 

[44] K. Glampedakis, S. A. Hughes, D. Kennefick, Phys. Rev. D 66, 044002 (2002), [arXiv:gr-qc/0203086]. 
[45] P.C. Peters, J. Matthews, Phys. Rev. 131, 435 (1963); P. C. Peters, Phys. Rev. 136, B1224 (1964). 
[46] L. Barack, A. Ori, Phys. Rev. D 64, 124003 (2001), [arXiv:gr-qc/0107056]. 
[47] P. L. Chrzanowski, Phys. Rev. D 11, 2042 (1975). 



[48] S. Hopper, C. R. Evans, Phys. Rev. D 82, 084010 (2010), [arXiv.org: 1006. 

[49] C. Cutler, D. Kennefick, E. Poisson, Phys. Rev. D 50, 3816 (1994). 

[50] A. Pound, Phys. Rev. D 81, 024023, (2010), [arXiv:0907.5107 [gr-qc]]. 

[51] E. Rosenthal, Phys. Rev. D 74, 084018 (2006). 

[52] [http://www.gnu.org/software/gsl/manual/html_node/index.html]. 

[53] L. Barack, private communication, 2010. 

[54] S. A. Teukolsky, Astrophys. J. 185, 635 (1973). 

[55] S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974). 



